Context Navigation

source: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx@ 352

Visit:

Last change on this file since 352 was 344, checked in by petry, 25 years ago
The pixelization in previous versions was buggy. This is the first version with a correct pixelization.
File size: 84.0 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.4 $
24	// $Author: petry $
25	// $Date: 2000-01-25 08:36:23 $
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
708	TFile outfile ( rootname , "RECREATE" ) ;
709
710	//
711	// create a Tree for the Event data stream
712	//
713
714	TTree EvtTree("EvtTree","Events of Run");
715
716	bsize=128000; split=1;
717
718	EvtTree.Branch("MRawEvt","MRawEvt",
719	&Evt, bsize, split);
720
721	EvtTree.Branch("MMcEvt","MMcEvt",
722	&McEvt, bsize, split);
723
724
725	float flli = 0. ;
726	unsigned short ulli = 0 ;
727
728	#endif // __ROOT__
729
730	// for safety and for dimensioning image_data: count the elements in the
731	// diagnostic data branch
732
733	i=0;
734	i++; // "n"
735	i++; // "primary"
736	i++; // "energy"
737	i++; // "cored"
738	i++; // "impact"
739	i++; // "xcore"
740	i++; // "ycore"
741	i++; // "theta"
742	i++; // "phi"
743	i++; // "deviations"
744	i++; // "dtheta"
745	i++; // "dphi"
746	i++; // "trigger"
747	i++; // "ncphs"
748	i++; // "maxpassthr_phe"
749	i++; // "nphes"
750	i++; // "nphes2"
751	i++; // "length"
752	i++; // "width"
753	i++; // "dist"
754	i++; // "xdist"
755	i++; // "azw"
756	i++; // "miss"
757	i++; // "alpha"
758	i++; // "conc2"
759	i++; // "conc3"
760	i++; // "conc4"
761	i++; // "conc5"
762	i++; // "conc6"
763	i++; // "conc7"
764	i++; // "conc8"
765	i++; // "conc9"
766	i++; // "conc10"
767	i++; // "asymx"
768	i++; // "asymy"
769	i++; // "phiasym"
770
771	nvar = i;
772	image_data = new float[nvar];
773
774	// set plate scale (deg/cm) and trigger area (deg)
775
776	plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;
777
778	if ( ! get_trigger_radius( &degTriggerZone ) )
779	degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
780
781	if ( ! get_correction( &fCorrection ) )
782	fCorrection = 1.0;
783
784	// number of pixels for parameters
785
786	anaPixels = get_ana_pixels();
787	anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
788
789	// open input file if we DO read data from a file
790
791	if (! Data_From_STDIN) {
792	log( SIGNATURE, "Openning input \"rfl\" file %s\n", inname );
793	inputfile.open( inname );
794	if ( inputfile.bad() )
795	error( SIGNATURE, "Cannot open input file: %s\n", inname );
796	}
797
798	// get signature, and check it
799	// NOTE: this part repeats further down in the code;
800	// if you change something here you probably want to change it
801	// there as well
802
803	strcpy(Signature, REFL_SIGNATURE);
804
805	strcpy(sign, Signature);
806
807	if ( Data_From_STDIN )
808	cin.read( (char *)sign, strlen(Signature));
809	else
810	inputfile.read( (char *)sign, strlen(Signature));
811
812	if (strcmp(sign, Signature) != 0) {
813	cerr << "ERROR: Signature of .rfl file is not correct\n";
814	cerr << '"' << sign << '"' << '\n';
815	cerr << "should be: " << Signature << '\n';
816	exit(1);
817	}
818
819	if ( Data_From_STDIN )
820	cin.read( (char *)sign, 1);
821	else
822	inputfile.read( (char *)sign, 1);
823
824	// open output file
825
826	log( SIGNATURE, "Openning output \"phe\" file %s\n", outname );
827	outputfile.open( outname );
828
829	if ( outputfile.bad() )
830	error( SIGNATURE, "Cannot open output file: %s\n", outname );
831
832	// open data file
833
834	log( SIGNATURE, "Openning data \"dat\" file %s\n", datname );
835	datafile.open( datname );
836
837	if ( outputfile.bad() )
838	error( SIGNATURE, "Cannot open output file: %s\n", outname );
839
840	// write signature
841
842	outputfile.write( SIGNATURE, sizeof(SIGNATURE) );
843
844	// initializes flag
845
846	strcpy( flag, " \0" );
847
848	// allocate space for PMTs numbers of pixels
849
850	fnpix = new float [ ct_NPixels ];
851	fnpixclean = new float [ ct_NPixels ];
852
853	#ifdef __ROOT__
854
855	fnslicesum = new float [ (2 * SLICES) ] ;
856
857	float slices [CAMERA_PIXELS][ (2 * SLICES) ] ;
858	float slices2 [CAMERA_PIXELS][ SLICES ] ;
859
860	float trans [ SLICES ] ;
861	#endif // __ROOT__
862
863
864	moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
865
866	//!@' @#### Main loop.
867	//@'
868
869	//begin my version
870
871	// get flag
872
873	if ( Data_From_STDIN )
874	cin.read( flag, SIZE_OF_FLAGS );
875	else
876	inputfile.read ( flag, SIZE_OF_FLAGS );
877
878	// loop over the file
879
880	still_in_loop = TRUE;
881
882	while (
883	((! Data_From_STDIN) && (! inputfile.eof()))
884	\|\|
885	(Data_From_STDIN && still_in_loop)
886	) {
887
888	// reading .rfl files
889	if(!isA( flag, FLAG_START_OF_RUN )){
890	error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
891	}
892	else { // found start of run
893	nshow=0;
894	// read next flag
895
896	if ( Data_From_STDIN )
897	cin.read( flag, SIZE_OF_FLAGS );
898	else
899	inputfile.read ( flag, SIZE_OF_FLAGS );
900
901	while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
902	/*!@'
903
904	For the case \|FLAG\_START\_OF\_EVENT\|,
905	we read each Cherenkov photon, and follow these steps:
906
907	@enumerate
908
909	@- Transform XY-coordinates to IJ-coordinates.
910
911	@- With this, we obtain the pixel where the photon hits.
912
913	@- Use the wavelength $\lambda$ and the table of QE, and
914	calculate the estimated (third order interpolated) quantum
915	efficiency for that photon. The photon can be rejected.
916
917	@- If accepted, then add to the pixel.
918
919	@endenumerate
920
921	In principle, we should stop here, and use another program to
922	'smear' the image, to add the Night Sky Background, and to
923	simulate the trigger logic, but we will make this program
924	quick and dirty, and include all here.
925
926	If we are reading PHE files, we jump to the point where the
927	pixelization process already has finished.
928
929	*/
930
931	++nshow;
932	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
933
934	// get MCEventHeader
935
936	if ( Data_From_STDIN )
937	cin.read( (char*)&mcevth, mcevth.mysize() );
938	else
939	mcevth.read( inputfile );
940
941	// calculate core distance and impact parameter
942
943	coreD = mcevth.get_core(&coreX, &coreY);
944
945	// calculate impact parameter (shortest distance betwee the original
946	// trajectory of the primary (assumed shower-axis) and the
947	// direction where the telescope points to
948	//
949	// we use the following equation, given that the shower core position
950	// is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
951	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
952	// and (l2,m2,n2)
953	//
954	// \| \|
955	// \| x1-x2 y1-y2 z1-z2 \|
956	// \| \|
957	// + \| l1 m1 n1 \|
958	// - \| \|
959	// \| l2 m2 n2 \|
960	// \| \|
961	// dist = ------------------------------------ ( > 0 )
962	// [ \|l1 m1\|2 \|m1 n1\|2 \|n1 l1\|2 ]1/2
963	// [ \| \| + \| \| + \| \| ]
964	// [ \|l2 m2\| \|m2 n2\| \|n2 l2\| ]
965	//
966	// playing a little bit, we get this reduced for in our case:
967	//
968	//
969	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
970	// [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
971	// 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
972
973	// read the direction of the incoming shower
974
975	thetashw = mcevth.get_theta();
976	phishw = mcevth.get_phi();
977
978	// calculate vector for shower
979
980	l1 = sin(thetashw)*cos(phishw);
981	m1 = sin(thetashw)*sin(phishw);
982	n1 = cos(thetashw);
983
984	// read the deviation of the telescope with respect to the shower
985
986	mcevth.get_deviations ( &thetaCT, &phiCT );
987
988	if ( (thetaCT == 0.) && (phiCT == 0.) ) {
989
990	// CT was looking to the source (both lines are parallel)
991	// therefore, we calculate the impact parameter as the distance
992	// between the CT axis and the core position
993
994	impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
995
996	} else {
997
998	// the shower comes off-axis
999
1000	// obtain with this the final direction of the CT
1001
1002	thetaCT += thetashw;
1003	phiCT += phishw;
1004
1005	// calculate vector for telescope
1006
1007	l2 = sin(thetaCT)*cos(phiCT);
1008	m2 = sin(thetaCT)*sin(phiCT);
1009	n2 = cos(thetaCT);
1010
1011	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1012	den = (SQR(l1m2 - l2m1) +
1013	SQR(m1n2 - m2n1) +
1014	SQR(n1l2 - n2l1));
1015	den = sqrt(den);
1016
1017	impactD = fabs(num)/den;
1018
1019	// fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ",
1020	// thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2,
1021	// num, den);
1022
1023	}
1024
1025	// clear camera
1026
1027	for ( i=0; i<ct_NPixels; ++i ){
1028
1029	fnpix[i] = 0.0;
1030	#ifdef __ROOT__
1031	for ( ii=0; ii<(2 * SLICES); ii++ ) {
1032	slices [i][ii] = 0 ;
1033	}
1034	#endif // __ROOT__
1035	}
1036
1037	ntcph +=ncph;
1038	ncph = 0;
1039
1040	#ifdef __DETAIL_TRIGGER__
1041	//
1042	// clear Trigger
1043	//
1044
1045	Trigger.Reset() ;
1046	#endif // __DETAIL_TRIGGER__
1047
1048	//- - - - - - - - - - - - - - - - - - - - - - - - -
1049	// read photons and "map" them into the pixels
1050	//--------------------------------------------------
1051
1052	// initialize CPhoton
1053
1054	cphoton.fill(0., 0., 0., 0., 0., 0., 0., 0.);
1055
1056	// read the photons data
1057
1058	if ( Data_From_STDIN )
1059	cin.read( flag, SIZE_OF_FLAGS );
1060	else
1061	inputfile.read ( flag, SIZE_OF_FLAGS );
1062
1063	// loop over the photons
1064
1065	t_ini = -99999;
1066
1067	while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
1068
1069	memcpy( (char*)&cphoton, flag, SIZE_OF_FLAGS );
1070
1071	if ( Data_From_STDIN )
1072	cin.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
1073	else
1074	inputfile.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
1075
1076	// increase number of photons
1077
1078	ncph++;
1079
1080	t = cphoton.get_t() ;
1081
1082	if(t_ini == -99999){ // this is the first photon we read from this event
1083	t_ini = t; // memorize time
1084	}
1085
1086	// The photons don't come in chronological order!
1087	// Put the first photon at the center of the array by adding the constant SLICES
1088
1089	t_chan = (int) ((t - t_ini )/ WIDTH_TIMESLICE ) + SLICES ;
1090
1091	if (t_chan > (2 * SLICES)){
1092	log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n",
1093	t_chan, (2 * SLICES), (2 * SLICES));
1094	t_chan = (2 * SLICES);
1095	}
1096	else if(t_chan < 0){
1097	log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n",
1098	t_chan, 0, 0);
1099	t_chan = 0;
1100	}
1101
1102	// Pixelization
1103
1104	cx = cphoton.get_x();
1105	cy = cphoton.get_y();
1106
1107	// get wavelength
1108
1109	last_wl = wl;
1110	wl = cphoton.get_wl();
1111
1112	// check if photon has valid wavelength and is inside outermost camera radius
1113
1114	if( (wl > 800.0) \|\| (wl < 290.0) \|\|
1115	(sqrt(cxcx + cycy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)) ){
1116
1117	// read next CPhoton
1118	if ( Data_From_STDIN )
1119	cin.read( flag, SIZE_OF_FLAGS );
1120	else
1121	inputfile.read ( flag, SIZE_OF_FLAGS );
1122
1123	// go to beginning of loop, the photon is lost
1124	continue;
1125
1126	}
1127
1128	// cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
1129
1130	nPMT = -1;
1131
1132	for(i=0; i<ct_NPixels; i++){
1133	if( bpoint_is_in_pix( cx, cy, i, &cam) ){
1134	nPMT = i;
1135	break;
1136	}
1137	}
1138
1139	if(nPMT==-1){// the photon is in none of the pixels
1140
1141	// read next CPhoton
1142	if ( Data_From_STDIN )
1143	cin.read( flag, SIZE_OF_FLAGS );
1144	else
1145	inputfile.read ( flag, SIZE_OF_FLAGS );
1146
1147	// go to beginning of loop, the photon is lost
1148	continue;
1149	}
1150
1151
1152
1153	#ifdef __QE__
1154
1155	//!@' @#### QE simulation.
1156	//@'
1157
1158	//+++
1159	// QE simulation
1160	//---
1161
1162	// find data point to be used in Lagrange interpolation (-> k)
1163
1164	qeptr = (float **)QE[nPMT];
1165
1166	FindLagrange(qeptr,k,wl);
1167
1168	// if random > quantum efficiency, reject it
1169
1170	qe = Lagrange(qeptr,k,wl) / 100.0;
1171
1172	// fprintf(stdout, "%f\n", qe);
1173
1174	if ( RandomNumber > qe ) {
1175
1176	// read next CPhoton
1177	if ( Data_From_STDIN )
1178	cin.read( flag, SIZE_OF_FLAGS );
1179	else
1180	inputfile.read ( flag, SIZE_OF_FLAGS );
1181
1182	// go to beginning of loop
1183	continue;
1184
1185	}
1186
1187	#endif // __QE__
1188
1189	//+++
1190	// Cphoton is accepted
1191	//---
1192
1193	// increase the number of Cphs. in the PMT, i.e.,
1194	// increase in one unit the counter of the photons
1195	// stored in the pixel nPMT
1196
1197	fnpix[nPMT] += 1.0;
1198
1199	#ifdef __ROOT__
1200	fnslicesum[t_chan] += 1.0 ;
1201	slices[nPMT][t_chan] += 1.0 ;
1202	#endif // __ROOT__
1203
1204	#ifdef __DETAIL_TRIGGER__
1205	//
1206	// fill the Trigger class with this phe
1207	//
1208	//
1209
1210	Trigger.Fill( nPMT, ( t - t_ini ) ) ;
1211	#endif // __DETAIL_TRIGGER__
1212
1213	// read next CPhoton
1214
1215	if ( Data_From_STDIN )
1216	cin.read( flag, SIZE_OF_FLAGS );
1217	else
1218	inputfile.read ( flag, SIZE_OF_FLAGS );
1219
1220	} // end while, i.e. found end of event
1221
1222	log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1223	ncph, ntcph);
1224
1225	// show number of photons
1226
1227	//cout << ncph << " photons read . . . " << endl << flush;
1228
1229	// skip it ?
1230
1231	for ( i=0; i<nSkip; ++i ) {
1232	if (Skip[i] == (nshow+ntshow)) {
1233	i = -1;
1234	break;
1235	}
1236	}
1237
1238	// if after the previous loop, the exit value of i is -1
1239	// then the shower number is in the list of showers to be
1240	// skipped
1241
1242	if (i == -1) {
1243	log(SIGNATURE, "\t\tskipped!\n");
1244	continue;
1245	}
1246
1247	/*!@'
1248
1249	After reading all the Cherenkov photons for a given event,
1250	we have in the table of number of photons for each pixel
1251	only the 'raw' amount of Cherenkov photons @$n_p@$. Now, we
1252	should take this number as the mean value of the
1253	distribution of photons in that pixel @$p@$, following a
1254	Poisson distribution.
1255
1256	@[ n_p \equiv \mu_p @]
1257
1258	and with this number the amount of light coming from the
1259	shower is calculated @$\hat{n}_p@$.
1260
1261	Then, we calculate the amount of Night Sky Background we
1262	must introduce in that pixel @$p@$. We calculate this using
1263	again a Poisson distribution with mean @$\mu_\mathrm{NSB}@$
1264	(defined in the \|camera.h\| file). The value of
1265	@$\mu_\mathrm{NSB}@$ is obtained from measurements. With this
1266	value, the amount of photons @$\hat{n}_\mathrm{NSB}@$ coming
1267	from the Night Sky Background is calculated.
1268
1269	Finally, the amount of photons for that pixels is:
1270	@[ \hat{n}_p^\mathrm{final} = \hat{n}_p + \hat{n}_\mathrm{NSB} @]
1271
1272	*/
1273
1274	// after reading all the photons, our camera is filled
1275
1276	if ( Select_Energy ) {
1277	if (( mcevth.get_energy() < Select_Energy_le ) \|\|
1278	( mcevth.get_energy() > Select_Energy_ue )) {
1279	log(SIGNATURE, "select_energy: shower rejected.\n");
1280	continue;
1281	}
1282	}
1283
1284	#ifdef __NSB__
1285
1286	//!@' @#### NSB (Night Sky Background) simulation.
1287	//@'
1288
1289	//+++
1290	// NSB simulation
1291	//---
1292
1293	// add NSB "noise"
1294	// TO DO: make meanNSB an array and read the contents from a file!
1295
1296	if ( simulateNSB )
1297	for ( i=0; i<ct_NPixels; ++i )
1298	fnpix[i] += (float)ignpoi( meanNSB );
1299
1300	#endif // __NSB__
1301
1302	// if we should apply any kind of correction, do it here.
1303
1304	for ( i=0; i<ct_NPixels; ++i )
1305	fnpix[i] *= fCorrection;
1306
1307	#ifdef __DETAIL_TRIGGER__
1308	// Trigger.Print() ;
1309	cout << Trigger.Diskriminate() << endl << endl ;
1310	#endif // __DETAIL_TRIGGER__
1311
1312	#ifdef __ROOT__
1313
1314	//
1315	// Fill the header of this event
1316	//
1317
1318	Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ;
1319
1320	// now put out the data of interest
1321	//
1322	// 1. -> look for the first slice with signal
1323	//
1324
1325	for ( i=0; i<(2 * SLICES); i++ )
1326	if ( fnslicesum[i] > 0. )
1327	break ;
1328
1329	startchan = i ;
1330
1331	//
1332	// copy the slices out of the big array
1333	//
1334	// put the first slice with signal to slice 4
1335	//
1336
1337	for (i=0; i<ct_NPixels; i++ )
1338	for ( ii=(startchan-3); ii < (startchan+12); ii++ )
1339	slices2 [i][ii-startchan+3] = slices [i][ii] ;
1340
1341
1342	//
1343	// if a pixes has a signal put it to the MRawEvt
1344	//
1345
1346	for (i=0; i<ct_NPixels; i++ ) {
1347	if ( fnpix[i] > 0 ) {
1348
1349	for ( ii=0; ii < 15; ii++ ) {
1350	trans [ii] = slices2[i][ii] ;
1351	}
1352
1353	Evt->FillPixel ( (UShort_t) i , trans ) ;
1354
1355	}
1356	}
1357
1358	//
1359	//
1360	//
1361
1362	McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1363	mcevth.get_energy(),
1364	mcevth.get_theta(),
1365	mcevth.get_phi(),
1366	mcevth.get_core(),
1367	mcevth.get_coreX(),
1368	mcevth.get_coreY(),
1369	flli,
1370	ulli, ulli, ulli, ulli, ulli ) ;
1371
1372	//
1373	// write it out to the file outfile
1374	//
1375
1376	EvtTree.Fill() ;
1377
1378	// clear all
1379
1380	Evt->Clear() ;
1381	McEvt->Clear() ;
1382
1383	#endif // __ROOT__
1384
1385	//++++++++++++++++++++++++++++++++++++++++++++++++++
1386	// at this point we have a camera full of
1387	// ph.e.s
1388	// we should first apply the trigger condition,
1389	// and if there's trigger, then clean the image,
1390	// calculate the islands statistics and the
1391	// other parameters of the image (Hillas' parameters
1392	// and so on).
1393	//--------------------------------------------------
1394
1395	#ifdef __TRIGGER__
1396
1397	/*!@'
1398
1399	@#### Trigger logic simulation.
1400
1401	In the following block we look at the pixel contents, looking
1402	for pixels fulfilling the trigger condition. This condition,
1403	in this current version of the program, is the following:
1404
1405	@itemize
1406
1407	@- \|CT1\|: Two neighbour pixels with charge above the threshold
1408	@$q_0@$. For the old CT1 data, however, the trigger condition
1409	was 'any two pixels with charge above the threshold @$q_0@$'.
1410
1411	@- \|MAGIC\|: A 'closed-packet' of four neighbour pixels, each
1412	of them with charge above the threshold @$q_0@$.
1413
1414	@enditemize
1415
1416	In the following figure you can find a sort of description
1417	about the meanning of 'closed-packet'.
1418
1419	@F
1420
1421	\begin{figure}[htbp]
1422	\begin{center}
1423	\includegraphics{closepck.eps}
1424	\caption{Meanning of the expression ``{\it close-packet}''}
1425	\label{fig:closepacket}
1426	\end{center}
1427	\end{figure}
1428
1429	@F
1430
1431	*/
1432
1433	//++
1434	// TRIGGER Condition
1435	//--
1436
1437	//@ If the input parameter "threshold" is 0 we find the maximum
1438	//@ trigger threshold this event can pass
1439
1440	for(k=0; ( qThreshold == 0 ? (k <= iMAX_THRESHOLD_PHE) : (k == 1) ); k++){
1441
1442	// is there trigger?
1443
1444	noverq0 = 0;
1445	q0 = ( qThreshold == 0. ? (float) k : qThreshold );
1446	trigger = FALSE;
1447	mxgrp = 0;
1448	maxcharge = 0.0;
1449
1450	// Warning! NOT all the camera is able to give trigger
1451	// only up to 'degTrigger' degrees
1452
1453	for ( i=0 ; (i<ct_NCentralPixels) && (trigger==FALSE) ; ++i ) {
1454
1455	// calculate absolute maximum
1456
1457	maxcharge = MAX(fnpix[i],maxcharge);
1458
1459	// is this pixel above threshold ?
1460
1461	if ( fnpix[i] <= q0 )
1462	continue;
1463
1464	// it is: increment the number of pixels above threshold
1465
1466	++noverq0;
1467
1468	// if the trigger already fired, just count the pixels
1469	// above threshold
1470
1471	if ( trigger == TRUE )
1472	continue;
1473
1474	// is this pixel inside the trigger zone in the camera ?
1475
1476	if ( (sqrt(SQR(pixary[i][0]) +
1477	SQR(pixary[i][1]))*plateScale_cm2deg) > degTriggerZone)
1478	continue;
1479
1480	// 'ngrpq0' is the number of neighbours of pixel i with q > q0
1481
1482	ngrpq0 = 0;
1483
1484	// look at each pixel in the neighborhood, and count
1485	// those above threshold q0
1486
1487	for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1 ; ++j )
1488	if ( fnpix[pixneig[i][j]] > q0 )
1489	++ngrpq0;
1490
1491	// check whether we have trigger
1492
1493	if ( ct_Type == 0 ) {
1494
1495	//++ >>>>> CT1 <<<<<
1496
1497	#ifdef __CT1_NO_NEIGHBOURS__
1498
1499	if ( noverq0 > 1 )
1500	trigger = TRUE;
1501
1502	#else
1503
1504	if ( ngrpq0 > 0 )
1505	trigger = TRUE;
1506
1507	#endif
1508
1509	//-- >>>>> CT1 <<<<<
1510
1511	} else {
1512
1513	//++ >>>>> MAGIC <<<<<
1514
1515	// (at least 4 packed with at least q0 phes)
1516
1517	// there are 3 cases
1518	// 1. zero, one or two neighbours have enough charge: no trigger
1519	// 2. five or six neighbours with enough charge: trigger? sure!!
1520	// 3. three or four neighbours with q > q0 : we must look
1521	// for 'closeness'.
1522
1523	switch ( ngrpq0 ) {
1524
1525	case 0:
1526	case 1:
1527	case 2:
1528
1529	trigger = FALSE;
1530	break;
1531
1532	case 3:
1533	case 4:
1534
1535	// if reaches this line, it means three or four neighbours
1536	// around the central pixel
1537
1538	triggerBits = 1;
1539
1540	for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) {
1541
1542	if ( fnpix[pixneig[i][j]] > q0 ) {
1543
1544	if ( pixary[pixneig[i][j]][0] > pixary[i][0] ) {
1545
1546	if ( nint(pixary[pixneig[i][j]][1]*10.0) >
1547	nint(pixary[i][1]*10.0) )
1548	bit = 2;
1549	else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
1550	nint(pixary[i][1]*10.0) )
1551	bit = 6;
1552	else
1553	bit = 1;
1554
1555	} else {
1556
1557	if ( nint(pixary[pixneig[i][j]][1]*10.0) >
1558	nint(pixary[i][1]*10.0) )
1559	bit = 3;
1560	else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
1561	nint(pixary[i][1]*10.0) )
1562	bit = 5;
1563	else
1564	bit = 4;
1565
1566	}
1567
1568	triggerBits \|= (1<<bit);
1569
1570	}
1571
1572	}
1573
1574	if ( ngrpq0 == 3 ) { // 4-fold trigger
1575
1576	switch ( triggerBits ) {
1577
1578	case 0x0f: // 0 000111 1
1579	case 0x1d: // 0 001110 1
1580	case 0x39: // 0 011100 1
1581	case 0x71: // 0 111000 1
1582	case 0x63: // 0 110001 1
1583	case 0x47: // 0 100011 1
1584
1585	trigger = TRUE;
1586	break;
1587
1588	default:
1589
1590	trigger = FALSE;
1591
1592	}
1593
1594	} else { // 4-fold trigger
1595
1596	switch ( triggerBits ) {
1597
1598	case 0x1f: // 0 001111 1
1599	case 0x3d: // 0 011110 1
1600	case 0x79: // 0 111100 1
1601	case 0x73: // 0 111001 1
1602	case 0x67: // 0 110011 1
1603	case 0x4f: // 0 100111 1
1604
1605	trigger = TRUE;
1606	break;
1607
1608	default:
1609
1610	trigger = FALSE;
1611
1612	}
1613
1614	}
1615
1616	mxgrp = MAX(ngrpq0,mxgrp);
1617
1618	break;
1619
1620	case 5:
1621	case 6:
1622
1623	trigger = TRUE;
1624	break;
1625
1626	default:
1627
1628	trigger = FALSE;
1629	error( SIGNATURE, "Number of neighbours > 6 !!! Exiting.\n\n");
1630	break;
1631
1632	} // switch (ngrpq0)
1633
1634	} // ct_Type
1635
1636	} // for each pixel i
1637
1638	if ( trigger == FALSE ) {
1639	break;
1640	} // end if
1641
1642	} //end for each threshold
1643	maxtrigthr_phe = (float) k-1; // i.e. maxtrigthr_phe < 0. means that
1644	// the event doesn't even pass threshold 0.
1645	// maxtrigthr_phe >= 0 means, the event passes some threshold
1646	// or (in case the input parameter "threshold" was > 0), the event
1647	// passes the threshold given by the input parameter.
1648	if ( maxtrigthr_phe >= 0. ) {
1649	trigger = TRUE;
1650	}
1651
1652
1653	novq0 = noverq0;
1654
1655	if ( trigger == TRUE ) {
1656
1657	itrigger = i;
1658	++ntrigger;
1659
1660	memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels );
1661
1662	#ifdef __TAILCUT__
1663
1664	//!@' @#### Tail-cut condition.
1665	//@'
1666
1667	//++
1668	// tail-cut
1669	//--
1670
1671	// Tail-Cut = 0 : No Tail-Cut
1672	// Tail-Cut > 0 : Make Tail-Cut
1673	// Tail-Cut < 0 : Make Tail-Cut with t_0 = Sqrt[ maximum ]
1674
1675	if (qTailCut > 0.0) {
1676
1677	for ( i=0; i<ct_NPixels; ++i )
1678	if ( fnpix[i] < qTailCut )
1679	fnpixclean[i] = 0.0;
1680
1681	} else if (qTailCut < 0.0) {
1682
1683	maxcharge = sqrt(maxcharge);
1684	for ( i=0; i<ct_NPixels; ++i )
1685	if ( fnpix[i] < maxcharge )
1686	fnpixclean[i] = 0.0;
1687
1688	}
1689
1690	#endif // __TAILCUT__
1691
1692	#ifdef __ISLANDS__
1693
1694	//!@' @#### Islands algorithm.
1695	//@'
1696
1697	//++
1698	// islands counting, and cleanning
1699	//--
1700
1701	if ( countIslands )
1702	do_islands( ct_NPixels, fnpixclean, pixneig, npixneig,
1703	countIslands, nIslandsCut);
1704
1705	#endif // __ISLANDS__
1706
1707	#ifdef __MOMENTS__
1708
1709	//!@' @#### Calculation of parameters of the image.
1710	//@'
1711
1712	//++
1713	// moments calculation
1714	//--
1715
1716	// calculate moments and other things
1717
1718	moments_ptr = moments( anaPixels, fnpixclean, pixary,
1719	plateScale_cm2deg, 0 );
1720
1721	charge = moments_ptr->charge ;
1722	smax = moments_ptr->smax ;
1723	maxs = moments_ptr->maxs ;
1724	nmaxs = moments_ptr->nmaxs ;
1725	length = moments_ptr->length ;
1726	width = moments_ptr->width ;
1727	dist = moments_ptr->dist ;
1728	xdist = moments_ptr->xdist ;
1729	azw = moments_ptr->azw ;
1730	miss = moments_ptr->miss ;
1731	alpha = moments_ptr->alpha ;
1732	conc = moments_ptr->conc ;
1733	asymx = moments_ptr->asymx ;
1734	asymx = moments_ptr->asymx ;
1735	phiasym= moments_ptr->phi;
1736
1737	lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary,
1738	plateScale_cm2deg,
1739	ct_PixelWidth_corner_2_corner_half);
1740
1741
1742	// fill the diagnostic Tree
1743
1744	event = new MDiagEventobject();
1745
1746	i=0;
1747	image_data[i] = event->n = hidt/10; i++;
1748	image_data[i] = event->primary = mcevth.get_primary(); i++;
1749	image_data[i] = event->energy = mcevth.get_energy(); i++;
1750	image_data[i] = event->cored = coreD; i++;
1751	image_data[i] = event->impact = impactD; i++;
1752	image_data[i] = event->xcore = coreX; i++;
1753	image_data[i] = event->ycore = coreY; i++;
1754	image_data[i] = event->theta = mcevth.get_theta(); i++;
1755	image_data[i] = event->phi = mcevth.get_phi(); i++;
1756	image_data[i] = event->deviations = mcevth.get_deviations (&dtheta, &dphi); i++;
1757	image_data[i] = event->dtheta = dtheta; i++;
1758	image_data[i] = event->dphi = dphi; i++;
1759	image_data[i] = event->trigger = trigger; i++;
1760	image_data[i] = event->ncphs = ncph; i++;
1761	image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++;
1762	image_data[i] = event->nphes = charge; i++;
1763	image_data[i] = event->nphes2 = smax; i++;
1764	image_data[i] = event->length = length; i++;
1765	image_data[i] = event->width = width; i++;
1766	image_data[i] = event->dist = dist; i++;
1767	image_data[i] = event->xdist = xdist; i++;
1768	image_data[i] = event->azw = azw; i++;
1769	image_data[i] = event->miss = miss; i++;
1770	image_data[i] = event->alpha = alpha; i++;
1771	image_data[i] = event->conc2 = conc[0]; i++;
1772	image_data[i] = event->conc3 = conc[1]; i++;
1773	image_data[i] = event->conc4 = conc[2]; i++;
1774	image_data[i] = event->conc5 = conc[3]; i++;
1775	image_data[i] = event->conc6 = conc[4]; i++;
1776	image_data[i] = event->conc7 = conc[5]; i++;
1777	image_data[i] = event->conc8 = conc[6]; i++;
1778	image_data[i] = event->conc9 = conc[7]; i++;
1779	image_data[i] = event->conc10 = conc[8]; i++;
1780	image_data[i] = event->asymx = asymx; i++;
1781	image_data[i] = event->asymy = asymy; i++;
1782	image_data[i] = event->phiasym = phiasym; i++;
1783
1784	// there should be "nvar" variables
1785
1786	if ( i != nvar )
1787	error( SIGNATURE, "Wrong entry number for diagnostic data.\n" );
1788
1789	tree->Fill();
1790	delete event;
1791
1792	// put information in the data file,
1793
1794	datafile << ntrigger;
1795	for(i=0;i<nvar;i++) {
1796	datafile << ' ' << image_data[i];
1797	}
1798
1799
1800	#endif // __MOMENTS__
1801
1802
1803	// revert the fnpixclean matrix into fnpix
1804	// (now we do this, but maybe in a future we want to
1805	// use both fnpix and fnpixclean for different things
1806
1807	memcpy( fnpix, fnpixclean, sizeof(float) * ct_NPixels );
1808
1809	// put this information in the data file,
1810
1811	if ( Write_All_Data ) {
1812	datafile << ' ' << -9999;
1813	for ( i=0; i<ct_NPixels; ++i )
1814	datafile << ' ' << fnpix[i];
1815	}
1816
1817	datafile << endl;
1818
1819	mcevth.set_trigger( TRUE );
1820
1821	log(SIGNATURE, "TRIGGER\n");
1822
1823	} else { // ( trigger == FALSE )
1824
1825	event = new MDiagEventobject();
1826
1827	i=0;
1828	image_data[i] = event->n = hidt/10; i++;
1829	image_data[i] = event->primary = mcevth.get_primary(); i++;
1830	image_data[i] = event->energy = mcevth.get_energy(); i++;
1831	image_data[i] = event->cored = coreD = mcevth.get_core(&coreX, &coreY); i++;
1832	image_data[i] = event->impact = coreD; i++;
1833	image_data[i] = event->xcore = coreX; i++;
1834	image_data[i] = event->ycore = coreY; i++;
1835	image_data[i] = event->theta = mcevth.get_theta(); i++;
1836	image_data[i] = event->phi = mcevth.get_phi(); i++;
1837	image_data[i] = event->deviations = mcevth.get_deviations(&dtheta, &dphi); i++;
1838	image_data[i] = event->dtheta = dtheta; i++;
1839	image_data[i] = event->dphi = dphi; i++;
1840	image_data[i] = event->trigger = trigger; i++;
1841	image_data[i] = event->ncphs = ncph; i++;
1842	image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++;
1843	image_data[i] = -1.; i++;
1844	image_data[i] = -1.; i++;
1845	image_data[i] = -1.; i++;
1846	image_data[i] = -1.; i++;
1847	image_data[i] = -1.; i++;
1848	image_data[i] = -1.; i++;
1849	image_data[i] = -1.; i++;
1850	image_data[i] = -1.; i++;
1851	image_data[i] = -1.; i++;
1852	image_data[i] = -1.; i++;
1853	image_data[i] = -1.; i++;
1854	image_data[i] = -1.; i++;
1855	image_data[i] = -1.; i++;
1856	image_data[i] = -1.; i++;
1857	image_data[i] = -1.; i++;
1858	image_data[i] = -1.; i++;
1859	image_data[i] = -1.; i++;
1860	image_data[i] = -1.; i++;
1861	image_data[i] = -1.; i++;
1862	image_data[i] = -1.; i++;
1863	image_data[i] = -1.; i++;
1864
1865	// there should be "nvar" variables
1866
1867	if ( i != nvar )
1868	error( SIGNATURE, "Wrong entry length for Ntuple.\n" );
1869
1870	tree->Fill();
1871	delete event;
1872
1873	// put this information in the data file,
1874
1875	if ( Write_All_Data ) {
1876
1877	datafile << ntrigger;
1878	for ( i=0; i<nvar; ++i )
1879	datafile << ' ' << image_data[i];
1880
1881	datafile << -9999;
1882	for ( i=0; i<ct_NPixels; ++i )
1883	datafile << ' ' << fnpix[i];
1884
1885	datafile << endl;
1886	}
1887
1888	mcevth.set_trigger( FALSE );
1889
1890	} // trigger == FALSE
1891
1892	#endif // __TRIGGER__
1893
1894	//!@' @#### Save data.
1895	//@'
1896
1897	//++++++++++++++++++++++++++++++++++++++++++++++++++
1898	// we now have all information we want
1899	// the only thing we must do now is writing it to
1900	// the output file
1901	//--------------------------------------------------
1902
1903	//++
1904	// save the image to the file
1905	//--
1906
1907	// write MCEventHeader to output file
1908
1909	outputfile.write( (char *)&mcevth, mcevth.mysize() );
1910
1911	#ifdef __TRIGGER__
1912
1913	// save the image
1914
1915	if ( (trigger == TRUE) \|\| (Write_All_Images == TRUE) )
1916	outputfile.write( (char ) fnpix, ct_NPixels sizeof( float ) );
1917
1918	#else
1919
1920	// save the image
1921
1922	outputfile.write( (char ) fnpix, ct_NPixels sizeof( float ) );
1923
1924	#endif // __TRIGGER__
1925
1926	if ( Data_From_STDIN )
1927	cin.read( flag, SIZE_OF_FLAGS );
1928	else
1929	inputfile.read ( flag, SIZE_OF_FLAGS );
1930
1931	} // end while there is a next event
1932
1933	if( !isA( flag, FLAG_END_OF_RUN )){
1934	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1935	}
1936	else { // found end of run
1937	ntshow += nshow;
1938	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1939
1940	if ( Data_From_STDIN )
1941	cin.read( flag, SIZE_OF_FLAGS );
1942	else
1943	inputfile.read ( flag, SIZE_OF_FLAGS );
1944
1945	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1946	log(SIGNATURE, "End of file . . .\n");
1947	still_in_loop = FALSE;
1948
1949	if ((! Data_From_STDIN) && (! inputfile.eof())){
1950
1951	// we have concatenated input files.
1952	// get signature of the next part and check it.
1953	// NOTE: this part repeats further up in the code;
1954	// if you change something here you probably want to change it
1955	// there as well
1956
1957	strcpy(Signature, REFL_SIGNATURE);
1958
1959	strcpy(sign, Signature);
1960
1961	inputfile.read( (char *)sign, strlen(Signature));
1962
1963	if (strcmp(sign, Signature) != 0) {
1964	cerr << "ERROR: Signature of .rfl file is not correct\n";
1965	cerr << '"' << sign << '"' << '\n';
1966	cerr << "should be: " << Signature << '\n';
1967	exit(1);
1968	}
1969
1970	if ( Data_From_STDIN )
1971	cin.read( (char *)sign, 1);
1972	else
1973	inputfile.read( (char *)sign, 1);
1974
1975	}
1976
1977	} // end if found end of file
1978	} // end if found end of run
1979	if ( Data_From_STDIN )
1980	cin.read( flag, SIZE_OF_FLAGS );
1981	else
1982	inputfile.read ( flag, SIZE_OF_FLAGS );
1983	} // end if else found start of run
1984	} // end big while loop
1985
1986	//!@' @#### End of program.
1987	//@'
1988
1989	//end my version
1990
1991	#ifdef __ROOT__
1992	//++
1993	// put the Event to the root file
1994	//--
1995
1996	EvtTree.Write() ;
1997	outfile.Write() ;
1998	outfile.Close() ;
1999
2000	#endif // __ROOT__
2001
2002	// close input file
2003
2004	ntcph += ncph;
2005	log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
2006	ntshow, ntcph );
2007	log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
2008	((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
2009
2010	// close files
2011
2012	log( SIGNATURE, "Closing files\n" );
2013
2014	inputfile.close();
2015	outputfile.close();
2016	datafile.close();
2017
2018	hfile->Write();
2019
2020	hfile->Close();
2021
2022	#ifdef __DETAIL_TRIGGER__
2023	// Output of Trigger statistics
2024	//
2025
2026	Trigger.PrintStat() ;
2027	#endif // __DETAIL_TRIGGER__
2028
2029	// program finished
2030
2031	log( SIGNATURE, "Done.\n");
2032
2033	return( 0 );
2034	}
2035	//!@}
2036
2037	// @T \newpage
2038
2039	//!@subsection Functions definition.
2040
2041	//!-----------------------------------------------------------
2042	// @name present
2043	//
2044	// @desc Make some presentation
2045	//
2046	// @date Sat Jun 27 05:58:56 MET DST 1998
2047	//------------------------------------------------------------
2048	// @function
2049
2050	//!@{
2051	void
2052	present(void)
2053	{
2054	cout << "##################################################\n"
2055	<< SIGNATURE << '\n' << '\n'
2056	<< "Processor of the reflector output\n"
2057	<< "J C Gonzalez, Jun 1998\n"
2058	<< "##################################################\n\n"
2059	<< flush ;
2060	}
2061	//!@}
2062
2063
2064	//!-----------------------------------------------------------
2065	// @name usage
2066	//
2067	// @desc show help
2068	//
2069	// @date Tue Dec 15 16:23:30 MET 1998
2070	//------------------------------------------------------------
2071	// @function
2072
2073	//!@{
2074	void
2075	usage(void)
2076	{
2077	present();
2078	cout << "\nusage ::\n\n"
2079	<< "\t camera "
2080	<< " [ -@ paramfile ] "
2081	<< " [ -h ] "
2082	<< "\n\n or \n\n"
2083	<< "\t camera < paramfile"
2084	<< "\n\n";
2085	exit(0);
2086	}
2087	//!@}
2088
2089
2090	//!-----------------------------------------------------------
2091	// @name log
2092	//
2093	// @desc function to send log information
2094	//
2095	// @var funct Name of the caller function
2096	// @var fmt Format to be used (message)
2097	// @var ... Other information to be shown
2098	//
2099	// @date Sat Jun 27 05:58:56 MET DST 1998
2100	//------------------------------------------------------------
2101	// @function
2102
2103	//!@{
2104	void
2105	log(const char funct, char fmt, ...)
2106	{
2107	va_list args;
2108
2109	// Display the name of the function that called error
2110	printf("[%s]: ", funct);
2111
2112	// Display the remainder of the message
2113	va_start(args, fmt);
2114	vprintf(fmt, args);
2115	va_end(args);
2116	}
2117	//!@}
2118
2119
2120	//!-----------------------------------------------------------
2121	// @name error
2122	//
2123	// @desc function to send an error message, and abort the program
2124	//
2125	// @var funct Name of the caller function
2126	// @var fmt Format to be used (message)
2127	// @var ... Other information to be shown
2128	//
2129	// @date Sat Jun 27 05:58:56 MET DST 1998
2130	//------------------------------------------------------------
2131	// @function
2132
2133	//!@{
2134	void
2135	error(const char funct, char fmt, ...)
2136	{
2137	va_list args;
2138
2139	// Display the name of the function that called error
2140	fprintf(stderr, "ERROR in %s: ", funct);
2141
2142	// Display the remainder of the message
2143	va_start(args, fmt);
2144	vfprintf(stderr, fmt, args);
2145	va_end(args);
2146
2147	perror(funct);
2148
2149	exit(1);
2150	}
2151	//!@}
2152
2153
2154	//!-----------------------------------------------------------
2155	// @name isA
2156	//
2157	// @desc returns TRUE(FALSE), if the flag is(is not) the given
2158	//
2159	// @var s1 String to be searched
2160	// @var flag Flag to compare with string s1
2161	// @return TRUE: both strings match; FALSE: oth.
2162	//
2163	// @date Wed Jul 8 15:25:39 MET DST 1998
2164	//------------------------------------------------------------
2165	// @function
2166
2167	//!@{
2168	int
2169	isA( char * s1, const char * flag ) {
2170	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
2171	}
2172	//!@}
2173
2174
2175	//!-----------------------------------------------------------
2176	// @name read_ct_file
2177	//
2178	// @desc read CT definition file
2179	//
2180	// @date Sat Jun 27 05:58:56 MET DST 1998
2181	//------------------------------------------------------------
2182	// @function
2183
2184	//!@{
2185	void
2186	read_ct_file(void)
2187	{
2188	char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
2189	char token[ITEM_MAX_LENGTH]; //@< a single token
2190	int i, j; //@< dummy counters
2191
2192	log( "read_ct_file", "start.\n" );
2193
2194	ifstream ctin ( ct_filename );
2195
2196	if ( ctin.bad() )
2197	error( "read_ct_file",
2198	"Cannot open CT def. file: %s\n", ct_filename );
2199
2200	// loop till the "end" directive is reached
2201
2202	while (!ctin.eof()) {
2203
2204	// get line from stdin
2205
2206	ctin.getline(line, LINE_MAX_LENGTH);
2207
2208	// look for each item at the beginning of the line
2209
2210	for (i=0; i<=define_mirrors; i++)
2211	if (strstr(line, CT_ITEM_NAMES[i]) == line)
2212	break;
2213
2214	// if it is not a valid line, just ignore it
2215
2216	if (i == define_mirrors+1)
2217	continue;
2218
2219	// case block for each directive
2220
2221	switch ( i ) {
2222
2223	case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2224
2225	// get focal distance
2226
2227	sscanf(line, "%s %d", token, &ct_Type);
2228
2229	log( "read_ct_file", "<Type of Telescope>: %s\n",
2230	((ct_Type==0) ? "CT1" : "MAGIC") );
2231
2232	break;
2233
2234	case focal_distance: // <focal distance> [cm]
2235
2236	// get focal distance
2237
2238	sscanf(line, "%s %f", token, &ct_Focal_mean);
2239
2240	log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2241
2242	break;
2243
2244	case focal_std: // s(focal distance) [cm]
2245
2246	// get focal distance
2247
2248	sscanf(line, "%s %f", token, &ct_Focal_std);
2249
2250	log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2251
2252	break;
2253
2254	case point_spread: // <point spread> [cm]
2255
2256	// get point spread
2257
2258	sscanf(line, "%s %f", token, &ct_PSpread_mean);
2259
2260	log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2261
2262	break;
2263
2264	case point_std: // s(point spread) [cm]
2265
2266	// get point spread
2267
2268	sscanf(line, "%s %f", token, &ct_PSpread_std);
2269
2270	log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2271
2272	break;
2273
2274	case adjustment_dev: // s(adjustment_dev) [cm]
2275
2276	// get point spread
2277
2278	sscanf(line, "%s %f", token, &ct_Adjustment_std);
2279
2280	log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2281
2282	break;
2283
2284	case black_spot: // radius of the black spot in the center [cm]
2285
2286	// get black spot radius
2287
2288	sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2289
2290	log( "read_ct_file", "Radius of the black spots: %f cm\n",
2291	ct_BlackSpot_rad);
2292
2293	break;
2294
2295	case r_mirror: // radius of the mirrors [cm]
2296
2297	// get radius of mirror
2298
2299	sscanf(line, "%s %f", token, &ct_RMirror);
2300
2301	log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2302
2303	break;
2304
2305	case n_mirrors: // number of mirrors
2306
2307	// get the name of the output_file from the line
2308
2309	sscanf(line, "%s %d", token, &ct_NMirrors);
2310
2311	log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2312
2313	break;
2314
2315	case camera_width: // camera width [cm]
2316
2317	// get the name of the ct_file from the line
2318
2319	sscanf(line, "%s %f", token, &ct_CameraWidth);
2320
2321	log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2322
2323	break;
2324
2325	case n_pixels: // number of pixels
2326
2327	// get the name of the output_file from the line
2328
2329	sscanf(line, "%s %d", token, &ct_NPixels);
2330
2331	log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2332
2333	break;
2334
2335	case n_centralpixels: // number of central pixels
2336
2337	// get the name of the output_file from the line
2338
2339	sscanf(line, "%s %d", token, &ct_NCentralPixels);
2340
2341	log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
2342
2343	break;
2344
2345	case n_gappixels: // number of gap pixels
2346
2347	// get the name of the output_file from the line
2348
2349	sscanf(line, "%s %d", token, &ct_NGapPixels);
2350
2351	log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
2352
2353	break;
2354
2355	case pixel_width: // pixel width [cm]
2356
2357	// get the name of the ct_file from the line
2358
2359	sscanf(line, "%s %f", token, &ct_PixelWidth);
2360
2361	ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
2362	ct_PixelWidth_corner_2_corner_half =
2363	ct_PixelWidth_corner_2_corner * 0.50;
2364	ct_Apot = ct_PixelWidth / 2;
2365	ct_2Apot = ct_Apot * 2.0;
2366
2367	log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2368
2369	break;
2370
2371	case define_mirrors: // read table with the parameters of the mirrors
2372
2373	log( "read_ct_file", "Table of mirrors data:\n" );
2374
2375	// check whether the number of mirrors was already set
2376
2377	if ( ct_NMirrors == 0 )
2378	error( "read_ct_file", "NMirrors was not set.\n" );
2379
2380	// allocate memory for paths list
2381
2382	log( "read_ct_file", "Allocating memory for ct_data\n" );
2383
2384	ct_data = new float*[ct_NMirrors];
2385
2386	for (i=0; i<ct_NMirrors; i++)
2387	ct_data[i] = new float[CT_NDATA];
2388
2389	// read data
2390
2391	log( "read_ct_file", "Reading mirrors data...\n" );
2392
2393	for (i=0; i<ct_NMirrors; i++)
2394	for (j=0; j<CT_NDATA; j++)
2395	ctin >> ct_data[i][j];
2396
2397	break;
2398
2399	} // switch ( i )
2400
2401	} // end while
2402
2403	// end
2404
2405	log( "read_ct_file", "done.\n" );
2406
2407	return;
2408	}
2409	//!@}
2410
2411
2412	//!-----------------------------------------------------------
2413	// @name read_pixels
2414	//
2415	// @desc read pixels data
2416	//
2417	// @date Fri Mar 12 16:33:34 MET 1999
2418	//------------------------------------------------------------
2419	// @function
2420
2421	//!@{
2422	void
2423	read_pixels(struct camera *pcam)
2424	{
2425	ifstream qefile;
2426	char line[LINE_MAX_LENGTH];
2427	int n, i, j, k;
2428	float qe;
2429
2430	//------------------------------------------------------------
2431	// first, pixels' coordinates
2432
2433	pcam->inumpixels = ct_NPixels;
2434	pcam->inumcentralpixels = ct_NCentralPixels;
2435	pcam->inumgappixels = ct_NGapPixels;
2436	pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
2437	pcam->dpixdiameter_cm = ct_PixelWidth;
2438
2439	// initialize pixel numbers
2440
2441	for ( i=0; i<PIX_ARRAY_SIDE; ++i )
2442	for ( j=0; j<PIX_ARRAY_SIDE; ++j )
2443	pixels[i][j][PIXNUM] = -1;
2444
2445	pixary = new float* [2*ct_NCentralPixels];
2446	for ( i=0; i<2*ct_NCentralPixels; ++i )
2447	pixary[i] = new float[2];
2448
2449	pixneig = new int* [ct_NCentralPixels];
2450	for ( i=0; i<ct_NCentralPixels; ++i ) {
2451	pixneig[i] = new int[6];
2452	for ( j=0; j<6; ++j )
2453	pixneig[i][j] = -1;
2454	}
2455
2456	npixneig = new int[ct_NCentralPixels];
2457	for ( i=0; i<ct_NCentralPixels; ++i )
2458	npixneig[i] = 0;
2459
2460	// generate all coordinates
2461
2462	igen_pixel_coordinates(pcam);
2463
2464	// transfer coordinates to the working arrays for
2465	// the central pixels
2466
2467	for(k=0; k<ct_NCentralPixels; k++){
2468
2469	i = (int) pcam->di[k];
2470	j = (int) pcam->dj[k];
2471
2472	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k;
2473	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k];
2474	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k];
2475
2476	pixary[k][0] = pcam->dxc[k];
2477	pixary[k][1] = pcam->dyc[k];
2478	}
2479
2480	// calculate tables of neighbours
2481
2482	#ifdef __DEBUG__
2483	for ( n=0 ; n<ct_NPixels ; ++n ) {
2484	cout << "Para el pixel " << n << ": ";
2485	for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2486	if ( pixels_are_neig(n,i) == TRUE ) {
2487	pixneig[n][npixneig[n]] = i;
2488	pixneig[i][npixneig[i]] = n;
2489	cout << i << ' ';
2490	++npixneig[n];
2491	++npixneig[i];
2492	}
2493	}
2494	cout << endl << flush;
2495	}
2496	#else // ! __DEBUG__
2497	for ( n=0 ; n<ct_NCentralPixels ; ++n )
2498	for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2499	if ( pixels_are_neig(n,i) == TRUE ) {
2500	pixneig[n][npixneig[n]] = i;
2501	pixneig[i][npixneig[i]] = n;
2502	++npixneig[n];
2503	++npixneig[i];
2504	}
2505	#endif // ! __DEBUG__
2506
2507	#ifdef __DEBUG__
2508	for ( n=0 ; n<ct_NPixels ; ++n ) {
2509	cout << n << ':';
2510	for ( j=0; j<npixneig[n]; ++j)
2511	cout << ' ' << pixneig[n][j];
2512	cout << endl << flush;
2513	}
2514	#endif // __DEBUG__
2515
2516	//------------------------------------------------------------
2517	// second, pixels' QE
2518
2519	// try to open the file
2520
2521	log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE);
2522
2523	qefile.open( QE_FILE );
2524
2525	// if it is wrong or does not exist, exit
2526
2527	if ( qefile.bad() )
2528	error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2529
2530	// read file
2531
2532	log("read_pixels", "Reading data . . .\n");
2533
2534	i=-1;
2535
2536	while ( ! qefile.eof() ) {
2537
2538	// get line from the file
2539
2540	qefile.getline(line, LINE_MAX_LENGTH);
2541
2542	// skip if comment
2543
2544	if ( *line == '#' )
2545	continue;
2546
2547	// if it is the first valid value, it is the number of QE data points
2548
2549	if ( i < 0 ) {
2550
2551	// get the number of datapoints
2552
2553	sscanf(line, "%d", &pointsQE);
2554
2555	// allocate memory for the table of QEs
2556
2557	QE = new float ** [ct_NPixels];
2558
2559	for ( i=0; i<ct_NPixels; ++i ) {
2560	QE[i] = new float * [2];
2561	QE[i][0] = new float[pointsQE];
2562	QE[i][1] = new float[pointsQE];
2563	}
2564
2565	QElambda = new float [pointsQE];
2566
2567	for ( i=0; i<pointsQE; ++i ) {
2568	qefile.getline(line, LINE_MAX_LENGTH);
2569	sscanf(line, "%f", &QElambda[i]);
2570	}
2571
2572	i=0;
2573
2574	continue;
2575	}
2576
2577	// get the values (num-pixel, num-datapoint, QE-value)
2578
2579	sscanf(line, "%d %d %f", &i, &j, &qe);
2580
2581	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2582	((j-1) < pointsQE) && ((j-1) > -1) ) {
2583	QE[i-1][0][j-1] = QElambda[j-1];
2584	QE[i-1][1][j-1] = qe;
2585	}
2586
2587	}
2588
2589	// close file
2590
2591	qefile.close();
2592
2593	// end
2594
2595	log("read_pixels", "Done.\n");
2596
2597	}
2598	//!@}
2599
2600
2601	//!-----------------------------------------------------------
2602	// @name pixels_are_neig
2603	//
2604	// @desc check whether two pixels are neighbours
2605	//
2606	// @var pix1 Number of the first pixel
2607	// @var pix2 Number of the second pixel
2608	// @return TRUE: both pixels are neighbours; FALSE: oth.
2609	//
2610	// @date Wed Sep 9 17:58:37 MET DST 1998
2611	//------------------------------------------------------------
2612	// @function
2613
2614	//!@{
2615	int
2616	pixels_are_neig(int pix1, int pix2)
2617	{
2618	if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2619	SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2620	> ct_PixelWidth_corner_2_corner )
2621	return ( FALSE );
2622	else
2623	return ( TRUE );
2624	}
2625	//!@}
2626
2627	//!-----------------------------------------------------------
2628	// @name igen_pixel_coordinates
2629	//
2630	// @desc generate the pixel center coordinates
2631	//
2632	// @var *pcam structure camera containing all the
2633	// camera information
2634	// @return total number of pixels
2635	//
2636	// DP
2637	//
2638	// @date Thu Oct 14 10:41:03 CEST 1999
2639	//------------------------------------------------------------
2640	// @function
2641
2642	//!@{
2643	/****** igen_pixel_coordinates() *******************************/
2644
2645	int igen_pixel_coordinates(struct camera *pcam) {
2646	/* generate pixel coordinates, return value is number of pixels */
2647
2648	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2649	float fsegment_fract;
2650	double dtsize;
2651	double dhsize;
2652	double dpsize;
2653	double dxfirst_pix;
2654	double dyfirst_pix;
2655	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
2656	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
2657
2658
2659	double dstartx, dstarty; /* for the gap pixels and outer pixels */
2660	int j, nrow;
2661
2662	dpsize = pcam->dpixdiameter_cm;
2663	dtsize = dpsize * sqrt(3.) / 2.;
2664	dhsize = dpsize / 2.;
2665
2666	/* Loop over central pixels to generate co-ordinates */
2667
2668	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2669
2670	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2671
2672	pcam->dpixsizefactor[ipixno-1] = 1.;
2673
2674	in = 0;
2675
2676	i = 0;
2677	itot_inside_ring = 0;
2678	iring_no = 0;
2679
2680	/* Calculate the number of pixels out to and including the ring containing pixel number */
2681	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2682
2683	while (itot_inside_ring == 0){
2684
2685	iN = 3(i(i+1)) + 1;
2686
2687	if (ipixno <= iN){
2688	iring_no = i;
2689	itot_inside_ring = iN;
2690	}
2691
2692	i++;
2693	}
2694
2695
2696	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2697
2698	ipix_in_ring = 0;
2699	for (i = 0; i < iring_no; ++i){
2700
2701	ipix_in_ring = ipix_in_ring + 6;
2702	}
2703
2704	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2705	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2706	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2707
2708	isegment = 0;
2709	fsegment_fract = 0.;
2710	if (iring_no > 0) {
2711
2712	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2713
2714	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2715
2716	}
2717
2718	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2719	/* pixel width (flat to flat) dpsize. */
2720
2721	dxfirst_pix = dpsize*iring_no;
2722	dyfirst_pix = 0.;
2723
2724	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2725
2726	ddxseg1 = - dhsize*iring_no;
2727	ddyseg1 = dtsize*iring_no;
2728	ddxseg2 = -dpsize*iring_no;
2729	ddyseg2 = 0.;
2730	ddxseg3 = ddxseg1;
2731	ddyseg3 = -ddyseg1;
2732	ddxseg4 = -ddxseg1;
2733	ddyseg4 = -ddyseg1;
2734	ddxseg5 = -ddxseg2;
2735	ddyseg5 = 0.;
2736	ddxseg6 = -ddxseg1;
2737	ddyseg6 = ddyseg1;
2738
2739	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2740	/* anti-clockwise around the ring by adding the segment to segment vectors. */
2741
2742	switch (isegment) {
2743
2744	case 0:
2745
2746	pcam->dxc[ipixno-1] = 0.;
2747	pcam->dyc[ipixno-1] = 0.;
2748
2749	case 1:
2750	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2751	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2752
2753	break;
2754
2755	case 2:
2756
2757	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2758	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2759
2760	break;
2761
2762	case 3:
2763
2764	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
2765	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
2766
2767	break;
2768
2769	case 4:
2770
2771	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
2772	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
2773
2774	break;
2775
2776	case 5:
2777
2778	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
2779	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2780
2781	break;
2782
2783	case 6:
2784
2785	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2786	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2787
2788	break;
2789
2790	default:
2791
2792	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2793	return(0);
2794
2795	} /* end switch */
2796
2797	} /* end for */
2798
2799	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2800	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2801
2802	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2803
2804	j = pcam->inumcentralpixels;
2805
2806	for(i=0; i<pcam->inumgappixels; i=i+6){
2807	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2808	pcam->dyc[j + i ] = dstarty;
2809	pcam->dpixsizefactor[j + i] = 1.;
2810	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2811	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2812	pcam->dpixsizefactor[j + i + 1] = 1.;
2813	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2814	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2815	pcam->dpixsizefactor[j + i+ 2] = 1.;
2816	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2817	pcam->dyc[j + i + 3] = dstarty;
2818	pcam->dpixsizefactor[j + i+ 3] = 1.;
2819	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2820	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2821	pcam->dpixsizefactor[j + i+ 4] = 1.;
2822	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2823	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2824	pcam->dpixsizefactor[j + i + 5] = 1.;
2825	} /* end for */
2826	} /* end if */
2827
2828	/* generate positions of the outer pixels */
2829
2830	if( pcam->inumbigpixels > 0 ){
2831
2832	j = pcam->inumcentralpixels + pcam->inumgappixels;
2833
2834	for(i=0; i<pcam->inumbigpixels; i++){
2835	pcam->dpixsizefactor[j + i] = 2.;
2836	}
2837
2838	in = 0;
2839
2840	nrow = (int) ceil(dstartx / 2. / dpsize);
2841
2842	while(in < pcam->inumbigpixels){
2843
2844	pcam->dxc[j + in] = dstartx + dpsize;
2845	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2846	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2847	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2848	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2849	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2850	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2851	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2852	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2853	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2854	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2855	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2856	for(i=1; i<nrow; i++){
2857	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2858	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2859	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2860	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2861	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2862	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2863	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2864	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2865	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
2866	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
2867	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
2868	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
2869	}
2870	in = in + 6 * nrow;
2871	dstartx = dstartx + 2. * dpsize;
2872	nrow = nrow + 1;
2873	} /* end while */
2874
2875	} /* end if */
2876
2877	/* generate the ij coordinates */
2878
2879	for(i=0; i<pcam->inumpixels; i++){
2880	pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize;
2881	pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60;
2882
2883	// fprintf(stdout, "%d %f %f %f %f %f\n",
2884	// i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i],
2885	// pcam->dpixsizefactor[i]);
2886
2887	}
2888
2889	return(pcam->inumpixels);
2890
2891	}
2892	//!@}
2893
2894	//!-----------------------------------------------------------
2895	// @name bpoint_is_in_pix
2896	//
2897	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
2898	//
2899	// @var *pcam structure camera containing all the
2900	// camera information
2901	// @var dx, dy point coordinates in centimeters
2902	// @var ipixnum pixel number (starting at 0)
2903	// @return TRUE if the point is inside the pixel, FALSE otherwise
2904	//
2905	// DP
2906	//
2907	// @date Thu Oct 14 16:59:04 CEST 1999
2908	//------------------------------------------------------------
2909	// @function
2910
2911	//!@{
2912
2913	/****** bpoint_is_in_pix() *************************************/
2914
2915	int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
2916	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
2917	/* the pixel is assumed to be a "closed set" */
2918
2919	double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
2920	double c, xx, yy; /* auxiliary variable */
2921
2922	b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
2923	a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
2924	c = 1. - 1./sqrt(3.);
2925	if((ipixnum < 0)\|\|(ipixnum >= pcam->inumpixels)){
2926	fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
2927	fprintf(stderr, "Exiting.\n");
2928	exit(203);
2929	}
2930	xx = dx - pcam->dxc[ipixnum];
2931	yy = dy - pcam->dyc[ipixnum];
2932
2933	if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) \|\|
2934	((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
2935	return(TRUE); /* point is inside */
2936	}
2937	else{
2938	return(FALSE); /* point is outside */
2939	}
2940	}
2941
2942	//!@}
2943
2944	//------------------------------------------------------------
2945	// @name dist_r_P
2946	//
2947	// @desc distance straight line r - point P
2948	//
2949	// @date Sat Jun 27 05:58:56 MET DST 1998
2950	// @function @code
2951	//------------------------------------------------------------
2952	// dist_r_P
2953	//
2954	// distance straight line r - point P
2955	//------------------------------------------------------------
2956
2957	float
2958	dist_r_P(float a, float b, float c,
2959	float l, float m, float n,
2960	float x, float y, float z)
2961	{
2962	return (
2963	sqrt((SQR((a-x)m-(b-y)l) +
2964	SQR((b-y)n-(c-z)m) +
2965	SQR((c-z)l-(a-x)n))/
2966	(SQR(l)+SQR(m)+SQR(n))
2967	)
2968	);
2969	}
2970	// @endcode
2971
2972
2973	//=------------------------------------------------------------
2974	//!@subsection Log of this file.
2975
2976	//!@{
2977	//
2978	// $Log: not supported by cvs2svn $
2979	// Revision 1.3 2000/01/20 18:22:17 petry
2980	// Found little bug which makes camera crash if it finds a photon
2981	// of invalid wavelength. This bug is now fixed and the range
2982	// of valid wavelengths extended to 290 - 800 nm.
2983	// This is in preparation for the NSB simulation to come.
2984	// Dirk
2985	//
2986	// Revision 1.2 1999/11/19 08:40:42 harald
2987	// Now it is possible to compile the camera programm under osf1.
2988	//
2989	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
2990	// This the starting point for CVS controlled further developments of the
2991	// camera program. The program was originally written by Jose Carlos.
2992	// But here you can find a "rootified" version to the program. This means
2993	// that there is no hbook stuff in it now. Also the output of the
2994	// program changed to the MagicRawDataFormat.
2995	//
2996	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
2997	//
2998	// In the following you can see the README file of that version:
2999	//
3000	// ==================================================
3001	//
3002	// Fri Oct 22 1999 D.P.
3003	//
3004	// The MAGIC Monte Carlo System
3005	//
3006	// Camera Simulation Programme
3007	// ---------------------------
3008	//
3009	// 1) Description
3010	//
3011	// This version is the result of the fusion of H.K.'s
3012	// root_camera which is described below (section 2)
3013	// and another version by D.P. which had a few additional
3014	// useful features.
3015	//
3016	// The version compiles under Linux with ROOT 2.22 installed
3017	// (variable ROOTSYS has to be set).
3018	//
3019	// Compile as before simply using "make" in the root_camera
3020	// directory.
3021	//
3022	// All features of H.K.'s root_camera were retained.
3023	//
3024	// Additional features of this version are:
3025	//
3026	// a) HBOOK is no longer used and all references are removed.
3027	//
3028	// b) Instead of HBOOK, the user is given now the possibility of
3029	// having Diagnostic data in ROOT format as a complement
3030	// to the ROOT Raw data.
3031	//
3032	// This data is written to the file which is determined by
3033	// the new input parameter "diag_file" in the camera parameter
3034	// file.
3035	//
3036	// All source code file belonging to this part have filenames
3037	// starting with "MDiag".
3038	//
3039	// The user can read the output file using the following commands
3040	// in an interactive ROOT session:
3041	//
3042	// root [0] .L MDiag.so
3043	// root [1] new TFile("diag.root");
3044	// root [2] new TTreeViewer("T");
3045	//
3046	// This brings up a viewer from which all variables of the
3047	// TTree can be accessed and histogrammed. This example
3048	// assumes that you have named the file "diag.root", that
3049	// you are using ROOT version 2.22 or later and that you have
3050	// the shared object library "MDiag.so" which is produced
3051	// by the Makefile along with the executable "camera".
3052	//
3053	// ! The contents of the so-called diag file is not yet fixed.
3054	// ! At the moment it is what J.C.G. used to put into the HBOOK
3055	// ! ntuple. In future versions the moments calculation can be
3056	// ! removed and the parameter list be modified correspondingly.
3057	//
3058	// c) Now concatenated reflector files can be read. This is useful
3059	// if you have run the reflector with different parameters but
3060	// you want to continue the analysis with all reflector data
3061	// going into ONE ROOT outputfile.
3062	//
3063	// The previous camera version contained a bug which made reading
3064	// of two or more concatenated reflector files impossible.
3065	//
3066	// d) The reflector output format was changed. It is now version
3067	// 0.4 .
3068	// The change solely consists in a shortening of the flag
3069	// definition in the file
3070	//
3071	// include-MC/MCCphoton.hxx
3072	//
3073	// ! IF YOU WANT TO READ REFLECTOR FORMAT 0.3, you can easily
3074	// ! do so by recompiling camera with the previous version of
3075	// ! include-MC/MCCphoton.hxx.
3076	//
3077	// The change was necessary for saving space and better
3078	// debugging. From now on, this format can be frozen.
3079	//
3080	// ! For producing reflector output in the new format, you
3081	// ! of course have to recompile your reflector with the
3082	// ! new include-MC/MCCphoton.hxx .
3083	//
3084	// e) A first version of the pixelization with the larger
3085	// outer pixels is implemented. THIS IS NOT YET FULLY
3086	// TESTED, but first rough tests show that it works
3087	// at least to a good approximation.
3088	//
3089	// The present version implements the camera outline
3090	// with 18 "gap-pixels" and 595 pixels in total as
3091	// shown in
3092	//
3093	// http://sarastro.ifae.es/internal/home/hardware/camera/numbering.ps
3094	//
3095	// This change involved
3096	//
3097	// (i) The file pixels.dat is no longer needed. Instead
3098	// the coordinates are generated by the program itself
3099	// (takes maybe 1 second). In the file
3100	//
3101	// pixel-coords.txt
3102	//
3103	// in the same directory as this README, you find a list
3104	// of the coordinates generated by this new routine. It
3105	// has the format
3106	//
3107	// number i j x y size-factor
3108	//
3109	// where i and j are J.C.G.'s so called biaxis hexagonal
3110	// coordinates (for internal use) and x and y are the
3111	// coordinates of the pixel centers in the standard camera
3112	// coordinate system in units of centimeters. The value
3113	// of "size-factor" determines the linear size of the pixel
3114	// relative to the central pixels.
3115	//
3116	// (ii) The magic.def file has two additional parameters
3117	// which give the number of central pixels and the
3118	// number of gap pixels
3119	//
3120	// (iii) In camera.h and camera.cxx several changes were
3121	// necessary, among them the introduction of several
3122	// new functions
3123	//
3124	// The newly suggested outline with asymmetric Winston cones
3125	// will be implemented in a later version.
3126	//
3127	// f) phe files can no longer be read since this contradicts
3128	// our philosophy that the analysis should be done with other
3129	// programs like e.g. EVITA and not with "camera" itself.
3130	// This possibility was removed.
3131	//
3132	// g) ROOT is no longer invoked with an interactive interface.
3133	// In this way, camera can better be run as a batch program and
3134	// it uses less memory.
3135	//
3136	// h) small changes concerning the variable "t_chan" were necessary in
3137	// order to avoid segmentation faults: The variable is used as an
3138	// index and it went sometimes outside the limits when camera
3139	// was reading proton data. This is because the reflector files
3140	// don't contain the photons in a chronological order and also
3141	// the timespread can be considerably longer that the foreseen
3142	// digitisation timespan. Please see the source code of camera.cxx
3143	// round about line 1090.
3144	//
3145	// j) several unused variables were removed, a few warning messages
3146	// occur when you compile camera.cxx but these can be ignored at
3147	// the moment.
3148	//
3149	// In general the program is of course not finished. It still needs
3150	// debugging, proper trigger simulation, simulation of the asymmetric
3151	// version of the outer pixels, proper NSB simulation, adaption of
3152	// the diag "ntuple" contents to our need and others small improvements.
3153	//
3154	// In the directory rfl-files there is now a file in reflector format 0.4
3155	// containing a single event produced by the starfiled adder. It has
3156	// a duration of 30 ns and represents the region around the Crab Nebula.
3157	// Using the enclosed input parameter file, camera should process this
3158	// file without problems.
3159	//
3160	// 2) The README for the previous version of root_camera
3161	//
3162	// README for a preliminary version of the
3163	// root_camera program.
3164	//
3165	// root_camera is based on the program "camera"of Jose Carlos
3166	// Gonzalez. It was changed in the way that only the pixelisation
3167	// and the distibution of the phe to the FADCs works in a
3168	// first version.
3169	//
3170	// Using the #undef command most possibilities of the orignal
3171	// program are switched of.
3172	//
3173	// The new parts are signed by
3174	//
3175	// - ROOT or __ROOT__
3176	// nearly all important codelines for ROOT output are enclosed
3177	// in structures like
3178	// #ifdef __ROOT__
3179	//
3180	// code
3181	//
3182	// #endif __ROOT__
3183	//
3184	// In same case the new lines are signed by a comment with the word
3185	// ROOT in it.
3186	//
3187	// For timing of the pulse some variable names are changed.
3188	// (t0, t1, t --> t_ini, t_fin, t_1st, t_chan,...)
3189	// Look also for this changes.
3190	//
3191	// For the new root-file is also a change in readparm-files
3192	//
3193	//
3194	// - __DETAIL_TRIGGER__
3195	//
3196	// This is for the implementation of the current work on trigger
3197	// studies. Because the class MTrigger is not well documented it
3198	// isn´t a part of this tar file. Only a dummy File exists.
3199	//
3200	//
3201	//
3202	// With all files in the archive, the root_camera program should run.
3203	//
3204	// A reflector file is in the directory rfl-files
3205	//
3206	// ==================================================
3207	//
3208	// From now on, use CVS for development!!!!
3209	//
3210	//
3211	//
3212	// Revision 1.3 1999/10/22 15:01:28 petry
3213	// version sent to H.K. and N.M. on Fri Oct 22 1999
3214	//
3215	// Revision 1.2 1999/10/22 09:44:23 petry
3216	// first synthesized version which compiles and runs without crashing;
3217	//
3218	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
3219	// first synthesised version
3220	//
3221	// Revision 1.13 1999/03/15 14:59:05 gonzalez
3222	// camera-1_1
3223	//
3224	// Revision 1.12 1999/03/02 09:56:10 gonzalez
3225	// * empty log message *
3226	//
3227	//
3228	//!@}
3229
3230	//=EOF

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

Download in other formats: