Context Navigation

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

Visit:

Last change on this file since 384 was 383, checked in by blanch, 25 years ago
A first FADC simulation and a trigger simulation are already implemented. The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software. A loop over trigger threshold and some corretcions in the time range where it looks for a trigger will be implemented as soon as possible.
File size: 71.8 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.7 $
24	// $Author: blanch $
25	// $Date: 2000-03-24 18:10:46 $
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 "TApplication.h"
52
53	#include "TFile.h"
54	#include "TTree.h"
55	#include "TBranch.h"
56	#include "TCanvas.h"
57
58	#include "MTrigger.hxx"
59	#include "MFadc.hxx"
60
61	#include "MRawEvt.h"
62	#include "MMcEvt.h"
63	#include "MMcTrig.hxx"
64
65	/*!@"
66
67	All the defines are located in the file \|camera.h\|.
68
69	@"*/
70
71	#include "camera.h"
72
73	//!@}
74
75	/*!@"
76
77	The following set of flags are used in time of compilation. They do
78	not affect directly the behaviour of the program at run-time
79	(though, of course, if you disconnected the option for
80	implementation of the Trigger logic, you will not be able to use any
81	trigger at all. The 'default' values mean default in the sense of
82	what you got from the server when you obtained this program.
83
84	@"*/
85
86	//!@{
87
88	// flag for debugging (default: OFF )
89	#define __DEBUG__
90	#undef __DEBUG__
91
92
93	//!@}
94
95	//=-----------------------------------------------------------
96	//!@subsection Definition of global variables.
97
98	/*!@"
99
100	Now we define some global variables with data about the telescope,
101	such as "focal distance", number of pixels/mirrors,
102	"size of the camera", and so on.
103
104	@"*/
105
106	/*!@"
107
108	Depending on the telescope we are using (CT1 or MAGIC), the
109	information stored in the definition file is different.
110	The variable \|ct_Type\| has the value 0 when we use
111	CT1, and 1 when we use MAGIC.
112
113	@"*/
114
115	//!@{
116	static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC
117	//!@}
118
119	/*!@"
120
121	And this is the information about the whole telescope.
122
123	@"*/
124
125	//!@{
126
127	// parameters of the CT (from the CT definition file)
128
129	////@: Focal distances [cm]
130	//static float *ct_Focal;
131
132	//@: Mean Focal distances [cm]
133	static float ct_Focal_mean;
134
135	//@: STDev. Focal distances [cm]
136	static float ct_Focal_std;
137
138	//@: Mean Point Spread function [cm]
139	static float ct_PSpread_mean;
140
141	//@: STDev. Point Spread function [cm]
142	static float ct_PSpread_std;
143
144	//@: STDev. Adjustmente deviation [cm]
145	static float ct_Adjustment_std;
146
147	//@: Radius of the Black Spot in mirror [cm]
148	static float ct_BlackSpot_rad;
149
150	//@: Radius of one mirror [cm]
151	static float ct_RMirror;
152
153	//@: Camera width [cm]
154	static float ct_CameraWidth;
155
156	//@: Pixel width [cm]
157	static float ct_PixelWidth;
158
159	//@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)
160	static float ct_PixelWidth_corner_2_corner;
161
162	//@: ct_PixelWidth_corner_2_corner / 2
163	static float ct_PixelWidth_corner_2_corner_half;
164
165	//@: Number of mirrors
166	static int ct_NMirrors = 0;
167
168	//@: Number of pixels
169	static int ct_NPixels;
170
171	//@: Number of pixels
172	static int ct_NCentralPixels;
173
174	//@: Number of pixels
175	static int ct_NGapPixels;
176
177	//@: ct_Apot = ct_PixelWidth / 2
178	static float ct_Apot;
179
180	//@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth
181	static float ct_2Apot;
182
183	//@: name of the CT definition file to use
184	static char ct_filename[256];
185
186	//@: list of showers to be skipped
187	static int *Skip;
188
189	//@: number of showers to be skipped
190	static int nSkip=0;
191
192	//@: flag: TRUE: data come from STDIN; FALSE: from file
193	static int Data_From_STDIN = FALSE;
194
195	//@: flag: TRUE: write all images to output; FALSE: only triggered showers
196	static int Write_All_Images = FALSE;
197
198	//@: flag: TRUE: write all data to output; FALSE: only triggered showers
199	static int Write_All_Data = FALSE;
200
201	static int Write_McEvt = TRUE;
202	static int Write_McTrig = FALSE;
203	static int Write_RawEvt = FALSE;
204
205	//@: flag: TRUE: selection on the energy
206	static int Select_Energy = TRUE;
207
208	//@: Lower edge of the selected energy range (in GeV)
209	static float Select_Energy_le = 0.0;
210
211	//@: Upper edge of the selected energy range (in GeV)
212	static float Select_Energy_ue = 100000.0;
213
214	//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
215	static int FADC_Scan = FALSE;
216
217	//@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't
218	static int Trigger_Scan = FALSE;
219
220	//!@}
221
222	/*!@"
223
224	The following double-pointer is a 2-dimensional table with information
225	about each pixel. The routine read_pixels will generate
226	the information for filling it using igen_pixel_coordinates().
227
228	@"*/
229
230	//!@{
231	// Pointer to a tables/Arrays with information about the pixels
232	// and data stored on them with information about the pixels
233
234
235	//@: coordinates x,y for each pixel
236	static float **pixary;
237
238	//@: indexes of pixels neighbours of a given one
239	static int **pixneig;
240
241	//@: number of neighbours a pixel have
242	static int *npixneig;
243
244	//@: contents of the pixels (ph.e.)
245	static float *fnpix;
246
247	//@: contents of the pixels (ph.e.) after cleanning
248	static float *fnpixclean;
249
250
251	//!@}
252
253	/*!@"
254
255	The following double-pointer is a 2-dimensional table with the
256	Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
257	of the wavelength @$\lambda@$. The routine \|read_pixels()\| will read
258	also this information from the file \|qe.dat\|.
259
260	@"*/
261
262	//!@{
263	// Pointer to a table with QE, number of datapoints, and wavelengths
264
265	//@: table of QE
266	static float ***QE;
267
268	//@: number of datapoints for the QE curve
269	static int pointsQE;
270
271	//@: table of QE
272	static float *QElambda;
273	//!@}
274
275	/*!@"
276
277	The following double-pointer is a 2-dimensional table with information
278	about each mirror in the dish. The routine \|read_ct_file()\| will read
279	this information from the CT definition file.
280
281	@"*/
282
283	//!@{
284	// Pointer to a table with the following info.:
285
286	static float **ct_data;
287
288	/*
289	* TYPE=0 (CT1)
290	* i s rho theta x y z thetan phin xn yn zn
291	*
292	* i : number of the mirror
293	* s : arc length [cm]
294	* rho : polar rho of the position of the center of the mirror [cm]
295	* theta : polar angle of the position of the center of the mirror [cm]
296	* x : x coordinate of the center of the mirror [cm]
297	* y : y coordinate of the center of the mirror [cm]
298	* z : z coordinate of the center of the mirror [cm]
299	* thetan : polar theta angle of the direction where the mirror points to
300	* phin : polar phi angle of the direction where the mirror points to
301	* xn : xn coordinate of the normal vector in the center (normalized)
302	* yn : yn coordinate of the normal vector in the center (normalized)
303	* zn : zn coordinate of the normal vector in the center (normalized)
304	*
305	* TYPE=1 (MAGIC)
306	* i f sx sy x y z thetan phin
307	*
308	* i : number of the mirror
309	* f : focal distance of that mirror
310	* sx : curvilinear coordinate of mirror's center in X[cm]
311	* sy : curvilinear coordinate of mirror's center in X[cm]
312	* x : x coordinate of the center of the mirror [cm]
313	* y : y coordinate of the center of the mirror [cm]
314	* z : z coordinate of the center of the mirror [cm]
315	* thetan : polar theta angle of the direction where the mirror points to
316	* phin : polar phi angle of the direction where the mirror points to
317	* xn : xn coordinate of the normal vector in the center (normalized)
318	* yn : yn coordinate of the normal vector in the center (normalized)
319	* zn : zn coordinate of the normal vector in the center (normalized)
320	*/
321	//!@}
322
323	/*!@"
324
325	We define a table into where random numbers will be stored.
326	The routines used for random number generation are provided by
327	\|RANLIB\| (taken from NETLIB, \|www.netlib.org\|), and by
328	the routine \|double drand48(void)\| (prototype defined in
329	\|stdlib.h\|) through the macro \|RandomNumber\| defined in
330	\|camera.h\|.
331
332	@"*/
333
334
335	//!@}
336
337	extern char FileName[];
338
339
340	//=-----------------------------------------------------------
341	// @subsection Main program.
342
343	//!@{
344
345	//++++++++++++++++++++++++++++++++++++++++
346	// MAIN PROGRAM
347	//----------------------------------------
348
349	int main(int argc, char **argv)
350	{
351
352	//!@' @#### Definition of variables.
353	//@'
354
355	char inname[256]; //@< input file name
356	char starfieldname[256]; //@< starfield input file name
357
358	char datname[256]; //@< data (ASCII) output file name
359	char diagname[256]; //@< diagnistic output file (ROOT format)
360
361	char rootname[256] ; //@< ROOT file name
362
363	char parname[256]; //@< parameters file name
364
365	char sign[20]; //@< initialize sign
366
367	char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
368
369	FILE *inputfile; //@< stream for the input file
370	ofstream datafile; //@< stream for the data file
371
372	MCEventHeader mcevth; //@< Event Header class (MC)
373	MCCphoton cphoton; //@< Cherenkov Photon class (MC)
374
375	Photoelectron *photoe = NULL; //@< array of the photoelectrons of one event
376	int inumphe; //@< number of photoelectrons in an event
377
378	float arrtmin_ns; //@ arrival time of the first photoelectron
379	float arrtmax_ns; //@ arrival time of the last photoelectron
380
381	float thetaCT, phiCT; //@< parameters of a given shower
382	float thetashw, phishw; //@< parameters of a given shower
383	float coreD, coreX, coreY; //@< core position and distance
384	float impactD; //@< impact parameter
385	float l1, m1, n1; //@< auxiliary variables
386	float l2, m2, n2; //@< auxiliary variables
387	float num, den; //@< auxiliary variables
388
389	int nshow=0; //@< partial number of shower in a given run
390	int ntshow=0; //@< total number of showers
391	int ncph=0; //@< partial number of photons in a given run
392	int ntcph=0; //@< total number of photons
393
394	int i, j, k; //@< simple counters
395
396	int simulateNSB; //@< Will we simulate NSB?
397	float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel)
398	float diffnsb_phepns[iMAXNUMPIX]; //@< diffuse NSB values for each pixel derived
399	//@< from meanNSB
400	float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS]; //@< non-diffuse nsb
401	//@< photoelectron rates
402	float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
403	EXTWAVEBAND1,
404	EXTWAVEBAND2,
405	EXTWAVEBAND3,
406	EXTWAVEBAND4,
407	EXTWAVEBAND5
408	};
409	float baseline_mv[iMAXNUMPIX]; //@< The baseline (mV) caused by the NSB; to be subtracted
410	//@< in order to simulate the preamps' AC coupling
411
412	float qThreshold; //@< Threshold value
413	float qTailCut; //@< Tail Cut value
414	int nIslandsCut; //@< Islands Cut value
415	int countIslands; //@< Will we count the islands?
416	int anaPixels;
417
418	float fCorrection; //@< Factor to apply to pixel values (def. 1.)
419
420	int trigger; //@< trigger flag
421	int itrigger; //@< index of pixel fired
422	int ntrigger = 0; //@< number of triggers in the whole file
423
424	float plateScale_cm2deg; //@< plate scale (deg/cm)
425	float degTriggerZone; //@< trigger area in the camera (radius, in deg.)
426
427	float dtheta, dphi; //@< deviations of CT from shower axis
428
429	int still_in_loop = FALSE;
430
431	float *image_data;
432	int nvar, hidt;
433
434	struct camera cam; // structure holding the camera definition
435
436
437	//!@' @#### Definition of variables for \|getopt()\|.
438	//@'
439
440	int ch, errflg = 0; //@< used by getopt
441
442	/*!@'
443
444	@#### Beginning of the program.
445
446	We start with the main program. First we (could) make some
447	presentation, and follows the reading of the parameters file (now
448	from the \|stdin\|), the reading of the CT parameters file, and the
449	creation of the output file, where the processed data will be
450	stored.
451
452	*/
453
454	//++
455	// START
456	//--
457
458	// make unbuffered output
459
460	cout.setf ( ios::stdio );
461
462	// parse command line options (see reflector.h)
463
464	parname[0] = '\0';
465
466	optarg = NULL;
467	while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
468	switch (ch) {
469	case 'f':
470	strcpy(parname, optarg);
471	break;
472	case 'h':
473	usage();
474	break;
475	default :
476	errflg++;
477	}
478
479	// show help if error
480
481	if ( errflg>0 )
482	usage();
483
484	// make some sort of presentation
485
486	present();
487
488	// read parameters file
489
490	if ( strlen(parname) < 1 )
491	readparam(NULL);
492	else
493	readparam(parname);
494
495	// read data from file or from STDIN?
496
497	Data_From_STDIN = get_data_from_stdin();
498
499	// write all images, even those without trigger?
500
501	Write_All_Images = get_write_all_images();
502
503	// write all data (i.e., ph.e.s in pixels)
504
505	Write_All_Data = get_write_all_data();
506
507	Write_McEvt = get_write_McEvt() ;
508	Write_McTrig = get_write_McTrig() ;
509	Write_RawEvt = get_write_RawEvt() ;
510
511	FADC_Scan = get_FADC_Scan();
512	Trigger_Scan = get_Trigger_Scan();
513
514	// get filenames
515
516	strcpy( inname, get_input_filename() );
517	strcpy( starfieldname, get_starfield_filename() );
518	strcpy( datname, get_data_filename() );
519	strcpy( diagname, get_diag_filename() );
520	strcpy( rootname, get_root_filename() );
521	strcpy( ct_filename, get_ct_filename() );
522
523	// get different parameters of the simulation
524
525	qThreshold = get_threshold();
526	qTailCut = get_tail_cut();
527	simulateNSB = get_nsb( &meanNSB );
528	countIslands = get_islands_cut( &nIslandsCut );
529
530	// get selections on the parameters
531
532	Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
533
534	// log filenames information
535
536	log(SIGNATURE,
537	"%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",
538	"Filenames",
539	"In", inname,
540	"Stars", starfieldname,
541	"CT", ct_filename,
542	"Data", datname,
543	"Diag", diagname,
544	"ROOT", rootname
545	);
546
547
548	// log flags information
549
550	log(SIGNATURE,
551	"%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
552	"Flags",
553	"Data_From_STDIN", ONoff(Data_From_STDIN),
554	"Write_All_Images", ONoff(Write_All_Images),
555	"Write_All_Data", ONoff(Write_All_Data));
556
557	// log flags information
558
559	log(SIGNATURE,
560	"%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
561	"Root output",
562	"Write_McEvt", ONoff(Write_McEvt),
563	"Write_McTrig", ONoff(Write_McTrig),
564	"Write_RawEvt", ONoff(Write_RawEvt));
565
566	// log parameters information
567
568	log(SIGNATURE,
569	"%s:\n\t%20s: %f\n\t%20s: %f\n\t%20s: %f %s\n\t%20s: %f %s\n",
570	"Parameters",
571	"q0 (Threshold)", qThreshold,
572	"t0 (Tail-cut)", qTailCut,
573	"NSB (phes/pixel)", meanNSB, ONoff(simulateNSB),
574	"i0 (Islands-cut)", nIslandsCut, ONoff(countIslands));
575
576	// log selections
577
578	log(SIGNATURE,
579	"%s:\n\t%20s: %s (%f:%f)\n",
580	"Selections:",
581	"Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
582
583	// set all random numbers seeds
584
585	setall( get_seeds(0), get_seeds(1) );
586
587	// get list of showers to evt. skip
588
589	nSkip = get_nskip_showers();
590
591	if (nSkip > 0) {
592	Skip = new int[ nSkip ];
593	get_skip_showers( Skip );
594
595	log(SIGNATURE, "There are some showers to skip:\n");
596	for (i=0; i<nSkip; ++i)
597	log(SIGNATURE, "\tshower # %d\n", Skip[i]);
598	}
599
600	// read parameters from the ct.def file
601
602	read_ct_file();
603
604	// read camera setup
605
606	read_pixels(&cam);
607
608	// allocate memory for the photoelectrons
609
610	photoe = new Photoelectron[iMAXNUMPHE];
611
612	Int_t Lev0, Lev1, Lev2 ;
613
614	// initialise ROOT
615
616	TROOT simple("simple", "MAGIC Telescope Monte Carlo");
617
618	// initialise instance of Trigger and FADC classes
619
620	MTrigger Trigger ; //@< A instance of the Class MTrigger
621
622	MMcTrig *McTrig = new MMcTrig() ;
623
624	MFadc fadc ; //@< A instance of the Class MFadc
625
626
627	// Prepare the raw data output
628
629	MRawEvt *Evt = new MRawEvt() ;
630	MMcEvt *McEvt = new MMcEvt ();
631
632	// initalize the ROOT file
633	//
634
635	TFile outfile ( rootname , "RECREATE" );
636
637	// create a Tree for the Event data stream
638
639	TTree EvtTree("EvtTree","Events of Run");
640
641	Int_t bsize=128000; Int_t split=1;
642
643	EvtTree.Branch("MRawEvt","MRawEvt",
644	&Evt, bsize, split);
645
646	EvtTree.Branch("MMcEvt","MMcEvt",
647	&McEvt, bsize, split);
648
649	EvtTree.Branch("MMcTrig","MMcTrig",
650	&McTrig, bsize, split);
651
652	unsigned short ulli = 0 ;
653
654	TApplication theAppTrigger("App", &argc, argv);
655
656	if (gROOT->IsBatch()) {
657	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
658	// return 1;
659	}
660
661	if(FADC_Scan){
662	TApplication theAppFadc("App", &argc, argv);
663
664	if (gROOT->IsBatch()) {
665	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
666	// return 1;
667	}
668	}
669
670	// for safety and for dimensioning image_data: count the elements in the
671	// diagnostic data branch
672
673	i=0;
674	i++; // "n"
675	i++; // "primary"
676	i++; // "energy"
677	i++; // "cored"
678	i++; // "impact"
679	i++; // "xcore"
680	i++; // "ycore"
681	i++; // "theta"
682	i++; // "phi"
683	i++; // "deviations"
684	i++; // "dtheta"
685	i++; // "dphi"
686	i++; // "trigger"
687	i++; // "ncphs"
688	i++; // "maxpassthr_phe"
689	i++; // "nphes"
690	i++; // "nphes2"
691	i++; // "length"
692	i++; // "width"
693	i++; // "dist"
694	i++; // "xdist"
695	i++; // "azw"
696	i++; // "miss"
697	i++; // "alpha"
698	i++; // "conc2"
699	i++; // "conc3"
700	i++; // "conc4"
701	i++; // "conc5"
702	i++; // "conc6"
703	i++; // "conc7"
704	i++; // "conc8"
705	i++; // "conc9"
706	i++; // "conc10"
707	i++; // "asymx"
708	i++; // "asymy"
709	i++; // "phiasym"
710
711	nvar = i;
712	image_data = new float[nvar];
713
714	// set plate scale (deg/cm) and trigger area (deg)
715
716	plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;
717
718	if ( ! get_trigger_radius( &degTriggerZone ) )
719	degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
720
721	if ( ! get_correction( &fCorrection ) )
722	fCorrection = 1.0;
723
724	// number of pixels for parameters
725
726	anaPixels = get_ana_pixels();
727	anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
728
729	// prepare the NSB simulation
730
731	if( simulateNSB ){
732
733	//
734	// Calculate the non-diffuse NSB photoelectron rates
735	//
736
737	k = produce_nsbrates( starfieldname,
738	&cam,
739	photoe, // only a dummy here
740	nsbrate_phepns );
741	if (k != 0){
742	cout << "Error when reading starfield... \nExiting.\n";
743	exit(1);
744	}
745
746	// calculate diffuse rate correcting for the pixel size
747
748	for(i=0; i<cam.inumpixels; i++){
749	diffnsb_phepns[i] = meanNSB *
750	cam.dpixsizefactor[i] * cam.dpixsizefactor[i];
751	}
752
753	}
754
755	//
756	// Read the reflector file with the Cherenkov data
757	//
758
759	// select input file
760
761	if ( Data_From_STDIN ) {
762
763	inputfile = stdin;
764
765	}
766	else{
767
768	log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
769	inputfile = fopen( inname, "r" );
770	if ( inputfile == NULL )
771	error( SIGNATURE, "Cannot open input file: %s\n", inname );
772
773	}
774
775	// get signature, and check it
776
777	if(check_reflector_file( inputfile )==FALSE){
778	exit(1);
779	}
780
781	// open data file
782
783	log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
784	datafile.open( datname );
785
786	if ( datafile.bad() )
787	error( SIGNATURE, "Cannot open data file: %s\n", datname );
788
789	// initializes flag
790
791	strcpy( flag, " \0" );
792
793	// allocate space for PMTs numbers of pixels
794
795	fnpix = new float [ ct_NPixels ];
796	fnpixclean = new float [ ct_NPixels ];
797
798	// initialize baseline
799
800	for(i=0; i<cam.inumpixels; i++){
801	baseline_mv[i] = 0.;
802	}
803
804	//!@' @#### Main loop.
805	//@'
806
807	// get flag
808
809	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
810
811	// loop over the file
812
813	still_in_loop = TRUE;
814
815	while (
816	((! Data_From_STDIN) && ( !feof(inputfile) ))
817	\|\|
818	(Data_From_STDIN && still_in_loop)
819	) {
820
821	// reading .rfl files
822	if(!isA( flag, FLAG_START_OF_RUN )){
823	error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
824	}
825	else { // found start of run
826
827	nshow=0;
828
829	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
830
831	while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
832
833	//
834	// Clear Trigger and Fadc
835	//
836	Trigger.Reset() ;
837	fadc.Reset() ;
838
839	++nshow;
840	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
841
842	// get MCEventHeader
843
844	fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
845
846	// calculate core distance and impact parameter
847
848	coreD = mcevth.get_core(&coreX, &coreY);
849
850	// calculate impact parameter (shortest distance betwee the original
851	// trajectory of the primary (assumed shower-axis) and the
852	// direction where the telescope points to
853	//
854	// we use the following equation, given that the shower core position
855	// is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
856	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
857	// and (l2,m2,n2)
858	//
859	// \| \|
860	// \| x1-x2 y1-y2 z1-z2 \|
861	// \| \|
862	// + \| l1 m1 n1 \|
863	// - \| \|
864	// \| l2 m2 n2 \|
865	// \| \|
866	// dist = ------------------------------------ ( > 0 )
867	// [ \|l1 m1\|2 \|m1 n1\|2 \|n1 l1\|2 ]1/2
868	// [ \| \| + \| \| + \| \| ]
869	// [ \|l2 m2\| \|m2 n2\| \|n2 l2\| ]
870	//
871	// playing a little bit, we get this reduced for in our case:
872	//
873	//
874	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
875	// [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
876	// 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
877
878	// read the direction of the incoming shower
879
880	thetashw = mcevth.get_theta();
881	phishw = mcevth.get_phi();
882
883	// calculate vector for shower
884
885	l1 = sin(thetashw)*cos(phishw);
886	m1 = sin(thetashw)*sin(phishw);
887	n1 = cos(thetashw);
888
889	// read the deviation of the telescope with respect to the shower
890
891	mcevth.get_deviations ( &thetaCT, &phiCT );
892
893	if ( (thetaCT == 0.) && (phiCT == 0.) ) {
894
895	// CT was looking to the source (both lines are parallel)
896	// therefore, we calculate the impact parameter as the distance
897	// between the CT axis and the core position
898
899	impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
900
901	} else {
902
903	// the shower comes off-axis
904
905	// obtain with this the final direction of the CT
906
907	thetaCT += thetashw;
908	phiCT += phishw;
909
910	// calculate vector for telescope
911
912	l2 = sin(thetaCT)*cos(phiCT);
913	m2 = sin(thetaCT)*sin(phiCT);
914	n2 = cos(thetaCT);
915
916	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
917	den = (SQR(l1m2 - l2m1) +
918	SQR(m1n2 - m2n1) +
919	SQR(n1l2 - n2l1));
920	den = sqrt(den);
921
922	impactD = fabs(num)/den;
923
924	// fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ",
925	// thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2,
926	// num, den);
927
928	}
929
930	// read the photons and produce the photoelectrons
931
932	k = produce_phes( inputfile,
933	&cam,
934	WAVEBANDBOUND1,
935	WAVEBANDBOUND6,
936	photoe, // will be changed by the function!
937	&inumphe, // important for later: the size of photoe[]
938	fnpix, // will be changed by the function!
939	&ncph, // will be changed by the function!
940	&arrtmin_ns, // will be changed by the function!
941	&arrtmax_ns // will be changed by the function!
942	);
943
944	if( k != 0 ){ // non-zero returnvalue means error
945	cout << "Exiting.\n";
946	exit(1);
947	}
948
949	log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
950	ncph, ntcph);
951
952	ntcph += ncph;
953
954	// skip it ?
955
956	for ( i=0; i<nSkip; ++i ) {
957	if (Skip[i] == (nshow+ntshow)) {
958	i = -1;
959	break;
960	}
961	}
962
963	// if after the previous loop, the exit value of i is -1
964	// then the shower number is in the list of showers to be
965	// skipped
966
967	if (i == -1) {
968	log(SIGNATURE, "\t\tskipped!\n");
969	continue;
970	}
971
972	// energy cut
973
974	if ( Select_Energy ) {
975	if (( mcevth.get_energy() < Select_Energy_le ) \|\|
976	( mcevth.get_energy() > Select_Energy_ue )) {
977	log(SIGNATURE, "select_energy: shower rejected.\n");
978	continue;
979	}
980	}
981
982	// NSB simulation
983
984	if(simulateNSB){
985
986	k = produce_nsb_phes( &arrtmin_ns, // will be changed by the function!
987	&arrtmax_ns, // will be changed by the function!
988	thetaCT,
989	&cam,
990	nsbrate_phepns,
991	diffnsb_phepns,
992	ext,
993	fnpix, // will be changed by the function!
994	photoe, // will be changed by the function!
995	&inumphe, // important for later: the size of photoe[]
996	baseline_mv // will be generated by the function
997	);
998
999	if( k != 0 ){ // non-zero returnvalue means error
1000	cout << "Exiting.\n";
1001	exit(1);
1002	}
1003
1004	}// end if(simulateNSB) ...
1005
1006
1007	// cout << arrtmin_ns << " " << arrtmax_ns << "\n";
1008	// for(i=0; i<cam.inumpixels; i++){
1009	// cout << i << " " << baseline_mv[i] <<"\n";
1010	// }
1011
1012	cout << "Total number of phes: " << inumphe <<endl;
1013
1014	// TRIGGER HERE
1015
1016	//
1017	// Put values of the analog signal for
1018	// each pixel from photoe array
1019	//
1020
1021	for(i=0;i<inumphe;i++){
1022	Trigger.FillShow(photoe[i].ipixnum,float((photoe[i].iarrtime_ns-arrtmin_ns)));
1023	fadc.Fill( photoe[i].ipixnum,(photoe[i].iarrtime_ns-arrtmin_ns) , Trigger.FillShow(photoe[i].ipixnum,float((photoe[i].iarrtime_ns-arrtmin_ns))));
1024	}
1025
1026	//
1027	// now the noise of the electronic
1028	// (preamps, optical transmission,..) is introduced.
1029	// This is done inside the class MTrigger by the method ElecNoise.
1030	//
1031	Trigger.ElecNoise() ;
1032	fadc.ElecNoise() ;
1033
1034
1035	Trigger.Diskriminate() ;
1036
1037	//
1038	// look if in all the signals in the trigger signal branch
1039	// is a possible Trigger. Therefore we habe to diskriminate all
1040	// the simulated analog signals (Method Diskriminate in class
1041	// MTrigger). We look simultanously for the moments at which
1042	// there are more than TRIGGER_MULTI pixels above the
1043	// CHANNEL_THRESHOLD.
1044	//
1045
1046	McTrig->SetZeroLevel( Lev0 = (Short_t) Trigger.ZeroLevel() ) ;
1047
1048	Lev1 = Lev2 = 0 ;
1049
1050	//
1051	// Start the First Level Trigger simulation
1052	//
1053
1054	if ( Lev0 > 0 ) {
1055	McTrig->SetFirstLevel ( Lev1 = Trigger.FirstLevel() ) ;
1056	}
1057
1058	// for( i=0; i<inumphe; i++){
1059	// cout << "phe " << photoe[i].ipixnum << " " << photoe[i].iarrtime_ns << "\n";
1060	// }
1061
1062	//if ( trigger>0 ) {
1063	if (Lev1>0){
1064
1065	itrigger = i;
1066	++ntrigger;
1067
1068	memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels );
1069	}
1070	//
1071	// Fill the header of this event
1072	//
1073
1074	Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ;
1075
1076	//
1077	// fill the MMcEvt with all information
1078	//
1079
1080	McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1081	mcevth.get_energy(),
1082	mcevth.get_theta(),
1083	mcevth.get_phi(),
1084	mcevth.get_core(),
1085	mcevth.get_coreX(),
1086	mcevth.get_coreY(),
1087	impactD,
1088	ulli, ulli,
1089	(UShort_t) ncph,
1090	ulli,
1091	(UShort_t) ncph) ;
1092
1093	// We don not count phtons out of the camera.
1094
1095	//
1096	// write it out to the file outfile
1097	//
1098
1099	EvtTree.Fill() ;
1100
1101
1102	//
1103	// if a first level trigger occurred, then
1104	// 1. do some other stuff (not implemented)
1105	// 2. start the gui tool
1106
1107	if(FADC_Scan){
1108	if ( Lev0 > 0 ) {
1109	fadc.ShowSignal( McEvt, (Float_t) 60. ) ;
1110	}
1111	}
1112
1113	if(Trigger_Scan){
1114	if ( Lev0 > 0 ) {
1115	Trigger.ShowSignal(McEvt) ;
1116	}
1117	}
1118
1119	// clear all
1120	Evt->Clear() ;
1121	McEvt->Clear() ;
1122	McTrig->Clear() ;
1123
1124
1125	//++++++++++++++++++++++++++++++++++++++++++++++++++
1126	// at this point we have a camera full of
1127	// ph.e.s
1128	// we should first apply the trigger condition,
1129	// and if there's trigger, then clean the image,
1130	// calculate the islands statistics and the
1131	// other parameters of the image (Hillas' parameters
1132	// and so on).
1133	//--------------------------------------------------
1134
1135	#ifdef __DEBUG__
1136	printf("\n");
1137
1138	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1139
1140	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1141
1142	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1143
1144	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1145
1146	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1147	(int)pixels[ici][icj][PIXNUM],
1148	pixels[ici][icj][PIXX],
1149	pixels[ici][icj][PIXY],
1150	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1151
1152	}
1153	}
1154	}
1155	}
1156
1157	for (i=0; i<ct_NPixels; ++i) {
1158	printf("%d (%d): ", i, npixneig[i]);
1159	for (j=0; j<npixneig[i]; ++i)
1160	printf(" %d", pixneig[i][j]);
1161	printf("\n");
1162	}
1163
1164	#endif // __DEBUG__
1165
1166
1167	//!@' @#### Save data.
1168	//@'
1169
1170	//++++++++++++++++++++++++++++++++++++++++++++++++++
1171	// we now have all information we want
1172	// the only thing we must do now is writing it to
1173	// the output file
1174	//--------------------------------------------------
1175
1176	//++
1177	// save the image to the file
1178	//--
1179
1180
1181
1182	// look for the next event
1183
1184	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1185
1186	} // end while there is a next event
1187
1188	if( !isA( flag, FLAG_END_OF_RUN )){
1189	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1190	}
1191	else { // found end of run
1192	ntshow += nshow;
1193	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1194
1195	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1196
1197	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1198	log(SIGNATURE, "End of file . . .\n");
1199	still_in_loop = FALSE;
1200
1201	if ((! Data_From_STDIN) && ( !feof(inputfile) )){
1202
1203	// we have concatenated input files.
1204	// get signature of the next part and check it.
1205
1206	if(check_reflector_file( inputfile )==FALSE){
1207	exit(1);
1208	}
1209
1210	}
1211
1212	} // end if found end of file
1213	} // end if found end of run
1214
1215	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1216
1217	} // end if else found start of run
1218	} // end big while loop
1219
1220	//++
1221	// put the Event to the root file
1222	//--
1223
1224	EvtTree.Write() ;
1225	outfile.Write() ;
1226	outfile.Close() ;
1227
1228
1229	// close input file
1230
1231	log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
1232	ntshow, ntcph );
1233	log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
1234	((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
1235
1236	// close files
1237
1238	log( SIGNATURE, "Closing files\n" );
1239
1240	if( ! Data_From_STDIN ){
1241	fclose( inputfile );
1242	}
1243	datafile.close();
1244
1245	// program finished
1246
1247	log( SIGNATURE, "Done.\n");
1248
1249	return( 0 );
1250	}
1251	//!@}
1252
1253	// @T \newpage
1254
1255	//!@subsection Functions definition.
1256
1257	//!-----------------------------------------------------------
1258	// @name present
1259	//
1260	// @desc Make some presentation
1261	//
1262	// @date Sat Jun 27 05:58:56 MET DST 1998
1263	//------------------------------------------------------------
1264	// @function
1265
1266	//!@{
1267	void
1268	present(void)
1269	{
1270	cout << "##################################################\n"
1271	<< SIGNATURE << '\n' << '\n'
1272	<< "Processor of the reflector output\n"
1273	<< "J C Gonzalez, Jun 1998\n"
1274	<< "##################################################\n\n"
1275	<< flush ;
1276	}
1277	//!@}
1278
1279
1280	//!-----------------------------------------------------------
1281	// @name usage
1282	//
1283	// @desc show help
1284	//
1285	// @date Tue Dec 15 16:23:30 MET 1998
1286	//------------------------------------------------------------
1287	// @function
1288
1289	//!@{
1290	void
1291	usage(void)
1292	{
1293	present();
1294	cout << "\nusage ::\n\n"
1295	<< "\t camera "
1296	<< " [ -@ paramfile ] "
1297	<< " [ -h ] "
1298	<< "\n\n or \n\n"
1299	<< "\t camera < paramfile"
1300	<< "\n\n";
1301	exit(0);
1302	}
1303	//!@}
1304
1305
1306	//!-----------------------------------------------------------
1307	// @name log
1308	//
1309	// @desc function to send log information
1310	//
1311	// @var funct Name of the caller function
1312	// @var fmt Format to be used (message)
1313	// @var ... Other information to be shown
1314	//
1315	// @date Sat Jun 27 05:58:56 MET DST 1998
1316	//------------------------------------------------------------
1317	// @function
1318
1319	//!@{
1320	void
1321	log(const char funct, char fmt, ...)
1322	{
1323	va_list args;
1324
1325	// Display the name of the function that called error
1326	printf("[%s]: ", funct);
1327
1328	// Display the remainder of the message
1329	va_start(args, fmt);
1330	vprintf(fmt, args);
1331	va_end(args);
1332	}
1333	//!@}
1334
1335
1336	//!-----------------------------------------------------------
1337	// @name error
1338	//
1339	// @desc function to send an error message, and abort the program
1340	//
1341	// @var funct Name of the caller function
1342	// @var fmt Format to be used (message)
1343	// @var ... Other information to be shown
1344	//
1345	// @date Sat Jun 27 05:58:56 MET DST 1998
1346	//------------------------------------------------------------
1347	// @function
1348
1349	//!@{
1350	void
1351	error(const char funct, char fmt, ...)
1352	{
1353	va_list args;
1354
1355	// Display the name of the function that called error
1356	fprintf(stdout, "ERROR in %s: ", funct);
1357
1358	// Display the remainder of the message
1359	va_start(args, fmt);
1360	vfprintf(stdout, fmt, args);
1361	va_end(args);
1362
1363	perror(funct);
1364
1365	exit(1);
1366	}
1367	//!@}
1368
1369
1370	//!-----------------------------------------------------------
1371	// @name isA
1372	//
1373	// @desc returns TRUE(FALSE), if the flag is(is not) the given
1374	//
1375	// @var s1 String to be searched
1376	// @var flag Flag to compare with string s1
1377	// @return TRUE: both strings match; FALSE: oth.
1378	//
1379	// @date Wed Jul 8 15:25:39 MET DST 1998
1380	//------------------------------------------------------------
1381	// @function
1382
1383	//!@{
1384	int
1385	isA( char * s1, const char * flag ) {
1386	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
1387	}
1388	//!@}
1389
1390
1391	//!-----------------------------------------------------------
1392	// @name read_ct_file
1393	//
1394	// @desc read CT definition file
1395	//
1396	// @date Sat Jun 27 05:58:56 MET DST 1998
1397	//------------------------------------------------------------
1398	// @function
1399
1400	//!@{
1401	void
1402	read_ct_file(void)
1403	{
1404	char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
1405	char token[ITEM_MAX_LENGTH]; //@< a single token
1406	int i, j; //@< dummy counters
1407
1408	log( "read_ct_file", "start.\n" );
1409
1410	ifstream ctin ( ct_filename );
1411
1412	if ( ctin.bad() )
1413	error( "read_ct_file",
1414	"Cannot open CT def. file: %s\n", ct_filename );
1415
1416	// loop till the "end" directive is reached
1417
1418	while (!ctin.eof()) {
1419
1420	// get line from stdin
1421
1422	ctin.getline(line, LINE_MAX_LENGTH);
1423
1424	// look for each item at the beginning of the line
1425
1426	for (i=0; i<=define_mirrors; i++)
1427	if (strstr(line, CT_ITEM_NAMES[i]) == line)
1428	break;
1429
1430	// if it is not a valid line, just ignore it
1431
1432	if (i == define_mirrors+1)
1433	continue;
1434
1435	// case block for each directive
1436
1437	switch ( i ) {
1438
1439	case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
1440
1441	// get focal distance
1442
1443	sscanf(line, "%s %d", token, &ct_Type);
1444
1445	log( "read_ct_file", "<Type of Telescope>: %s\n",
1446	((ct_Type==0) ? "CT1" : "MAGIC") );
1447
1448	break;
1449
1450	case focal_distance: // <focal distance> [cm]
1451
1452	// get focal distance
1453
1454	sscanf(line, "%s %f", token, &ct_Focal_mean);
1455
1456	log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
1457
1458	break;
1459
1460	case focal_std: // s(focal distance) [cm]
1461
1462	// get focal distance
1463
1464	sscanf(line, "%s %f", token, &ct_Focal_std);
1465
1466	log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
1467
1468	break;
1469
1470	case point_spread: // <point spread> [cm]
1471
1472	// get point spread
1473
1474	sscanf(line, "%s %f", token, &ct_PSpread_mean);
1475
1476	log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
1477
1478	break;
1479
1480	case point_std: // s(point spread) [cm]
1481
1482	// get point spread
1483
1484	sscanf(line, "%s %f", token, &ct_PSpread_std);
1485
1486	log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
1487
1488	break;
1489
1490	case adjustment_dev: // s(adjustment_dev) [cm]
1491
1492	// get point spread
1493
1494	sscanf(line, "%s %f", token, &ct_Adjustment_std);
1495
1496	log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
1497
1498	break;
1499
1500	case black_spot: // radius of the black spot in the center [cm]
1501
1502	// get black spot radius
1503
1504	sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
1505
1506	log( "read_ct_file", "Radius of the black spots: %f cm\n",
1507	ct_BlackSpot_rad);
1508
1509	break;
1510
1511	case r_mirror: // radius of the mirrors [cm]
1512
1513	// get radius of mirror
1514
1515	sscanf(line, "%s %f", token, &ct_RMirror);
1516
1517	log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
1518
1519	break;
1520
1521	case n_mirrors: // number of mirrors
1522
1523	// get the name of the output_file from the line
1524
1525	sscanf(line, "%s %d", token, &ct_NMirrors);
1526
1527	log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
1528
1529	break;
1530
1531	case camera_width: // camera width [cm]
1532
1533	// get the name of the ct_file from the line
1534
1535	sscanf(line, "%s %f", token, &ct_CameraWidth);
1536
1537	log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
1538
1539	break;
1540
1541	case n_pixels: // number of pixels
1542
1543	// get the name of the output_file from the line
1544
1545	sscanf(line, "%s %d", token, &ct_NPixels);
1546
1547	log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
1548
1549	break;
1550
1551	case n_centralpixels: // number of central pixels
1552
1553	// get the name of the output_file from the line
1554
1555	sscanf(line, "%s %d", token, &ct_NCentralPixels);
1556
1557	log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
1558
1559	break;
1560
1561	case n_gappixels: // number of gap pixels
1562
1563	// get the name of the output_file from the line
1564
1565	sscanf(line, "%s %d", token, &ct_NGapPixels);
1566
1567	log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
1568
1569	break;
1570
1571	case pixel_width: // pixel width [cm]
1572
1573	// get the name of the ct_file from the line
1574
1575	sscanf(line, "%s %f", token, &ct_PixelWidth);
1576
1577	ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
1578	ct_PixelWidth_corner_2_corner_half =
1579	ct_PixelWidth_corner_2_corner * 0.50;
1580	ct_Apot = ct_PixelWidth / 2;
1581	ct_2Apot = ct_Apot * 2.0;
1582
1583	log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
1584
1585	break;
1586
1587	case define_mirrors: // read table with the parameters of the mirrors
1588
1589	log( "read_ct_file", "Table of mirrors data:\n" );
1590
1591	// check whether the number of mirrors was already set
1592
1593	if ( ct_NMirrors == 0 )
1594	error( "read_ct_file", "NMirrors was not set.\n" );
1595
1596	// allocate memory for paths list
1597
1598	log( "read_ct_file", "Allocating memory for ct_data\n" );
1599
1600	ct_data = new float*[ct_NMirrors];
1601
1602	for (i=0; i<ct_NMirrors; i++)
1603	ct_data[i] = new float[CT_NDATA];
1604
1605	// read data
1606
1607	log( "read_ct_file", "Reading mirrors data...\n" );
1608
1609	for (i=0; i<ct_NMirrors; i++)
1610	for (j=0; j<CT_NDATA; j++)
1611	ctin >> ct_data[i][j];
1612
1613	break;
1614
1615	} // switch ( i )
1616
1617	} // end while
1618
1619	// end
1620
1621	log( "read_ct_file", "done.\n" );
1622
1623	return;
1624	}
1625	//!@}
1626
1627
1628	//!-----------------------------------------------------------
1629	// @name read_pixels
1630	//
1631	// @desc read pixels data
1632	//
1633	// @date Fri Mar 12 16:33:34 MET 1999
1634	//------------------------------------------------------------
1635	// @function
1636
1637	//!@{
1638	void
1639	read_pixels(struct camera *pcam)
1640	{
1641	ifstream qefile;
1642	char line[LINE_MAX_LENGTH];
1643	int n, i, j, icount;
1644	float qe;
1645
1646	//------------------------------------------------------------
1647	// first, pixels' coordinates
1648
1649	pcam->inumpixels = ct_NPixels;
1650	pcam->inumcentralpixels = ct_NCentralPixels;
1651	pcam->inumgappixels = ct_NGapPixels;
1652	pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
1653	pcam->dpixdiameter_cm = ct_PixelWidth;
1654
1655	// initialize pixel numbers
1656
1657	pixary = new float* [2*ct_NCentralPixels];
1658	for ( i=0; i<2*ct_NCentralPixels; ++i )
1659	pixary[i] = new float[2];
1660
1661	pixneig = new int* [ct_NCentralPixels];
1662	for ( i=0; i<ct_NCentralPixels; ++i ) {
1663	pixneig[i] = new int[6];
1664	for ( j=0; j<6; ++j )
1665	pixneig[i][j] = -1;
1666	}
1667
1668	npixneig = new int[ct_NCentralPixels];
1669	for ( i=0; i<ct_NCentralPixels; ++i )
1670	npixneig[i] = 0;
1671
1672	// generate all coordinates
1673
1674	igen_pixel_coordinates(pcam);
1675
1676
1677	// calculate tables of neighbours
1678
1679	#ifdef __DEBUG__
1680	for ( n=0 ; n<ct_NPixels ; ++n ) {
1681	cout << "Para el pixel " << n << ": ";
1682	for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
1683	if ( pixels_are_neig(n,i) == TRUE ) {
1684	pixneig[n][npixneig[n]] = i;
1685	pixneig[i][npixneig[i]] = n;
1686	cout << i << ' ';
1687	++npixneig[n];
1688	++npixneig[i];
1689	}
1690	}
1691	cout << endl << flush;
1692	}
1693	#else // ! __DEBUG__
1694	for ( n=0 ; n<ct_NCentralPixels ; ++n )
1695	for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
1696	if ( pixels_are_neig(n,i) == TRUE ) {
1697	pixneig[n][npixneig[n]] = i;
1698	pixneig[i][npixneig[i]] = n;
1699	++npixneig[n];
1700	++npixneig[i];
1701	}
1702	#endif // ! __DEBUG__
1703
1704	#ifdef __DEBUG__
1705	for ( n=0 ; n<ct_NPixels ; ++n ) {
1706	cout << n << ':';
1707	for ( j=0; j<npixneig[n]; ++j)
1708	cout << ' ' << pixneig[n][j];
1709	cout << endl << flush;
1710	}
1711	#endif // __DEBUG__
1712
1713	//------------------------------------------------------------
1714	// second, pixels' QE
1715
1716	// try to open the file
1717
1718	log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
1719
1720	qefile.open( QE_FILE );
1721
1722	// if it is wrong or does not exist, exit
1723
1724	if ( qefile.bad() )
1725	error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
1726
1727	// read file
1728
1729	log("read_pixels", "Reading data . . .\n");
1730
1731	i=-1;
1732	icount = 0;
1733
1734	while ( ! qefile.eof() ) {
1735
1736	// get line from the file
1737
1738	qefile.getline(line, LINE_MAX_LENGTH);
1739
1740	// skip if comment
1741
1742	if ( *line == '#' )
1743	continue;
1744
1745	// if it is the first valid value, it is the number of QE data points
1746
1747	if ( i < 0 ) {
1748
1749	// get the number of datapoints
1750
1751	sscanf(line, "%d", &pointsQE);
1752
1753	// allocate memory for the table of QEs
1754
1755	QE = new float ** [ct_NPixels];
1756
1757	for ( i=0; i<ct_NPixels; ++i ) {
1758	QE[i] = new float * [2];
1759	QE[i][0] = new float[pointsQE];
1760	QE[i][1] = new float[pointsQE];
1761	}
1762
1763	QElambda = new float [pointsQE];
1764
1765	for ( i=0; i<pointsQE; ++i ) {
1766	qefile.getline(line, LINE_MAX_LENGTH);
1767	sscanf(line, "%f", &QElambda[i]);
1768	}
1769
1770	i=0;
1771
1772	continue;
1773	}
1774
1775	// get the values (num-pixel, num-datapoint, QE-value)
1776
1777	if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
1778	break;
1779
1780	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
1781	((j-1) < pointsQE) && ((j-1) > -1) ) {
1782	QE[i-1][0][j-1] = QElambda[j-1];
1783	QE[i-1][1][j-1] = qe;
1784	}
1785
1786	if ( i > ct_NPixels) break;
1787
1788	icount++;
1789
1790	}
1791
1792	if(icount/pointsQE < ct_NPixels){
1793	error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
1794	icount/pointsQE, ct_NPixels );
1795	}
1796
1797	// close file
1798
1799	qefile.close();
1800
1801	// test QE
1802
1803	for(icount=0; icount< ct_NPixels; icount++){
1804	for(i=0; i<pointsQE; i++){
1805	if( QE[icount][0][i] < 100. \|\| QE[icount][0][i] > 1000. \|\|
1806	QE[icount][1][i] < 0. \|\| QE[icount][1][i] > 100.){
1807	error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
1808	icount, i, QE[icount][0][i], QE[icount][1][i] );
1809	}
1810	}
1811	}
1812
1813	// end
1814
1815	log("read_pixels", "Done.\n");
1816
1817	}
1818	//!@}
1819
1820
1821	//!-----------------------------------------------------------
1822	// @name pixels_are_neig
1823	//
1824	// @desc check whether two pixels are neighbours
1825	//
1826	// @var pix1 Number of the first pixel
1827	// @var pix2 Number of the second pixel
1828	// @return TRUE: both pixels are neighbours; FALSE: oth.
1829	//
1830	// @date Wed Sep 9 17:58:37 MET DST 1998
1831	//------------------------------------------------------------
1832	// @function
1833
1834	//!@{
1835	int
1836	pixels_are_neig(int pix1, int pix2)
1837	{
1838	if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
1839	SQR( pixary[pix1][1] - pixary[pix2][1] ) )
1840	> ct_PixelWidth_corner_2_corner )
1841	return ( FALSE );
1842	else
1843	return ( TRUE );
1844	}
1845	//!@}
1846
1847	//!-----------------------------------------------------------
1848	// @name igen_pixel_coordinates
1849	//
1850	// @desc generate the pixel center coordinates
1851	//
1852	// @var *pcam structure camera containing all the
1853	// camera information
1854	// @return total number of pixels
1855	//
1856	// DP
1857	//
1858	// @date Thu Oct 14 10:41:03 CEST 1999
1859	//------------------------------------------------------------
1860	// @function
1861
1862	//!@{
1863	/****** igen_pixel_coordinates() *******************************/
1864
1865	int igen_pixel_coordinates(struct camera *pcam) {
1866	/* generate pixel coordinates, return value is number of pixels */
1867
1868	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
1869	float fsegment_fract;
1870	double dtsize;
1871	double dhsize;
1872	double dpsize;
1873	double dxfirst_pix;
1874	double dyfirst_pix;
1875	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
1876	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
1877
1878
1879	double dstartx, dstarty; /* for the gap pixels and outer pixels */
1880	int j, nrow;
1881
1882	dpsize = pcam->dpixdiameter_cm;
1883	dtsize = dpsize * sqrt(3.) / 2.;
1884	dhsize = dpsize / 2.;
1885
1886	/* Loop over central pixels to generate co-ordinates */
1887
1888	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
1889
1890	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
1891
1892	pcam->dpixsizefactor[ipixno-1] = 1.;
1893
1894	in = 0;
1895
1896	i = 0;
1897	itot_inside_ring = 0;
1898	iring_no = 0;
1899
1900	/* Calculate the number of pixels out to and including the ring containing pixel number */
1901	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
1902
1903	while (itot_inside_ring == 0){
1904
1905	iN = 3(i(i+1)) + 1;
1906
1907	if (ipixno <= iN){
1908	iring_no = i;
1909	itot_inside_ring = iN;
1910	}
1911
1912	i++;
1913	}
1914
1915
1916	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
1917
1918	ipix_in_ring = 0;
1919	for (i = 0; i < iring_no; ++i){
1920
1921	ipix_in_ring = ipix_in_ring + 6;
1922	}
1923
1924	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
1925	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
1926	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
1927
1928	isegment = 0;
1929	fsegment_fract = 0.;
1930	if (iring_no > 0) {
1931
1932	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
1933
1934	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
1935
1936	}
1937
1938	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
1939	/* pixel width (flat to flat) dpsize. */
1940
1941	dxfirst_pix = dpsize*iring_no;
1942	dyfirst_pix = 0.;
1943
1944	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
1945
1946	ddxseg1 = - dhsize*iring_no;
1947	ddyseg1 = dtsize*iring_no;
1948	ddxseg2 = -dpsize*iring_no;
1949	ddyseg2 = 0.;
1950	ddxseg3 = ddxseg1;
1951	ddyseg3 = -ddyseg1;
1952	ddxseg4 = -ddxseg1;
1953	ddyseg4 = -ddyseg1;
1954	ddxseg5 = -ddxseg2;
1955	ddyseg5 = 0.;
1956	ddxseg6 = -ddxseg1;
1957	ddyseg6 = ddyseg1;
1958
1959	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
1960	/* anti-clockwise around the ring by adding the segment to segment vectors. */
1961
1962	switch (isegment) {
1963
1964	case 0:
1965
1966	pcam->dxc[ipixno-1] = 0.;
1967	pcam->dyc[ipixno-1] = 0.;
1968
1969	case 1:
1970	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
1971	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
1972
1973	break;
1974
1975	case 2:
1976
1977	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
1978	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
1979
1980	break;
1981
1982	case 3:
1983
1984	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
1985	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
1986
1987	break;
1988
1989	case 4:
1990
1991	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
1992	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
1993
1994	break;
1995
1996	case 5:
1997
1998	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
1999	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2000
2001	break;
2002
2003	case 6:
2004
2005	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2006	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2007
2008	break;
2009
2010	default:
2011
2012	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2013	return(0);
2014
2015	} /* end switch */
2016
2017	} /* end for */
2018
2019	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2020	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2021
2022	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2023
2024	j = pcam->inumcentralpixels;
2025
2026	for(i=0; i<pcam->inumgappixels; i=i+6){
2027	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2028	pcam->dyc[j + i ] = dstarty;
2029	pcam->dpixsizefactor[j + i] = 1.;
2030	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2031	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2032	pcam->dpixsizefactor[j + i + 1] = 1.;
2033	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2034	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2035	pcam->dpixsizefactor[j + i+ 2] = 1.;
2036	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2037	pcam->dyc[j + i + 3] = dstarty;
2038	pcam->dpixsizefactor[j + i+ 3] = 1.;
2039	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2040	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2041	pcam->dpixsizefactor[j + i+ 4] = 1.;
2042	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2043	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2044	pcam->dpixsizefactor[j + i + 5] = 1.;
2045	} /* end for */
2046	} /* end if */
2047
2048	/* generate positions of the outer pixels */
2049
2050	if( pcam->inumbigpixels > 0 ){
2051
2052	j = pcam->inumcentralpixels + pcam->inumgappixels;
2053
2054	for(i=0; i<pcam->inumbigpixels; i++){
2055	pcam->dpixsizefactor[j + i] = 2.;
2056	}
2057
2058	in = 0;
2059
2060	nrow = (int) ceil(dstartx / 2. / dpsize);
2061
2062	while(in < pcam->inumbigpixels){
2063
2064	pcam->dxc[j + in] = dstartx + dpsize;
2065	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2066	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2067	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2068	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2069	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2070	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2071	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2072	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2073	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2074	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2075	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2076	for(i=1; i<nrow; i++){
2077	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2078	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2079	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2080	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2081	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2082	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2083	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2084	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2085	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
2086	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
2087	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
2088	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
2089	}
2090	in = in + 6 * nrow;
2091	dstartx = dstartx + 2. * dpsize;
2092	nrow = nrow + 1;
2093	} /* end while */
2094
2095	} /* end if */
2096
2097	return(pcam->inumpixels);
2098
2099	}
2100	//!@}
2101
2102	//!-----------------------------------------------------------
2103	// @name bpoint_is_in_pix
2104	//
2105	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
2106	//
2107	// @var *pcam structure camera containing all the
2108	// camera information
2109	// @var dx, dy point coordinates in centimeters
2110	// @var ipixnum pixel number (starting at 0)
2111	// @return TRUE if the point is inside the pixel, FALSE otherwise
2112	//
2113	// DP
2114	//
2115	// @date Thu Oct 14 16:59:04 CEST 1999
2116	//------------------------------------------------------------
2117	// @function
2118
2119	//!@{
2120
2121	/****** bpoint_is_in_pix() *************************************/
2122
2123	int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
2124	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
2125	/* the pixel is assumed to be a "closed set" */
2126
2127	double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
2128	double c, xx, yy; /* auxiliary variable */
2129
2130	b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
2131	a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
2132	c = 1. - 1./sqrt(3.);
2133	if((ipixnum < 0)\|\|(ipixnum >= pcam->inumpixels)){
2134	fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
2135	fprintf(stderr, "Exiting.\n");
2136	exit(203);
2137	}
2138	xx = dx - pcam->dxc[ipixnum];
2139	yy = dy - pcam->dyc[ipixnum];
2140
2141	if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) \|\|
2142	((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
2143	return(TRUE); /* point is inside */
2144	}
2145	else{
2146	return(FALSE); /* point is outside */
2147	}
2148	}
2149
2150	//!@}
2151
2152	//------------------------------------------------------------
2153	// @name dist_r_P
2154	//
2155	// @desc distance straight line r - point P
2156	//
2157	// @date Sat Jun 27 05:58:56 MET DST 1998
2158	// @function @code
2159	//------------------------------------------------------------
2160	// dist_r_P
2161	//
2162	// distance straight line r - point P
2163	//------------------------------------------------------------
2164
2165	float
2166	dist_r_P(float a, float b, float c,
2167	float l, float m, float n,
2168	float x, float y, float z)
2169	{
2170	return (
2171	sqrt((SQR((a-x)m-(b-y)l) +
2172	SQR((b-y)n-(c-z)m) +
2173	SQR((c-z)l-(a-x)n))/
2174	(SQR(l)+SQR(m)+SQR(n))
2175	)
2176	);
2177	}
2178
2179	//------------------------------------------------------------
2180	// @name check_reflector_file
2181	//
2182	// @desc check if a given reflector file has the right signature
2183	// @desc return TRUE or FALSE
2184	//
2185	// @date Mon Feb 14 16:44:21 CET 2000
2186	// @function @code
2187	//------------------------------------------------------------
2188
2189	int check_reflector_file(FILE *infile){
2190
2191	char Signature[20]; // auxiliary variable
2192	char sign[20]; // auxiliary variable
2193
2194	strcpy(Signature, REFL_SIGNATURE);
2195
2196	strcpy(sign, Signature);
2197
2198	fread( (char *)sign, strlen(Signature), 1, infile);
2199
2200	if (strcmp(sign, Signature) != 0) {
2201	cout << "ERROR: Signature of .rfl file is not correct\n";
2202	cout << '"' << sign << '"' << '\n';
2203	cout << "should be: " << Signature << '\n';
2204	return(FALSE);
2205	}
2206
2207	fread( (char *)sign, 1, 1, infile);
2208
2209	return(TRUE);
2210
2211	}
2212
2213	//------------------------------------------------------------
2214	// @name lin_interpol
2215	//
2216	// @desc interpolate linearly between two points returning the
2217	// @desc y-value of the result
2218	//
2219	// @date Thu Feb 17 11:31:32 CET 2000
2220	// @function @code
2221	//------------------------------------------------------------
2222
2223	float lin_interpol( float x1, float y1, float x2, float y2, float x){
2224
2225	if( (x2 - x1)<=0. ){ // avoid division by zero, return average
2226	cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
2227	return((y1+y2)/2.);
2228	}
2229	else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
2230	return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
2231	}
2232	}
2233
2234
2235	//------------------------------------------------------------
2236	// @name Photoelectron
2237	//
2238	// @desc constructor for class Photoelectron
2239	//
2240	// @date Mon Feb 15 16:44:21 CET 2000
2241	// @function @code
2242	//------------------------------------------------------------
2243
2244	Photoelectron::Photoelectron(void){
2245	iarrtime_ns = NOTIME;
2246	ipixnum = -1;
2247	}
2248
2249	//------------------------------------------------------------
2250	// @name produce_phes
2251	//
2252	// @desc read the photons of an event, pixelize them and simulate QE
2253	// @desc return various statistics and the array of Photoelectrons
2254	//
2255	// @date Mon Feb 14 16:44:21 CET 2000
2256	// @function @code
2257	//------------------------------------------------------------
2258
2259	int produce_phes( FILE *sp, // the input file
2260	struct camera *cam, // the camera layout
2261	float minwl_nm, // the minimum accepted wavelength
2262	float maxwl_nm, // the maximum accepted wavelength
2263	class Photoelectron phe[iMAXNUMPHE], // the generated phes
2264	int *itotnphe, // total number of produced photoelectrons
2265	float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
2266	int *incph, // total number of cph read
2267	float *tmin_ns, // minimum arrival time of all phes
2268	float *tmax_ns // maximum arrival time of all phes
2269	){
2270
2271	static int i, k, ipixnum;
2272	static float cx, cy, wl, qe, t;
2273	static MCCphoton photon;
2274	static float **qept;
2275	static char flag[SIZE_OF_FLAGS + 1];
2276	static float radius;
2277
2278	// reset variables
2279
2280	for ( i=0; i<cam->inumpixels; ++i ){
2281
2282	nphe[i] = 0.0;
2283
2284	}
2285
2286	*itotnphe = 0;
2287	*incph = 0;
2288	*tmin_ns = NOTIME; // very big
2289	*tmax_ns = -NOTIME; // very small
2290
2291	radius = cam->dxc[cam->inumpixels-1]
2292	+ 1.5cam->dpixdiameter_cmcam->dpixsizefactor[cam->inumpixels-1];
2293
2294	//- - - - - - - - - - - - - - - - - - - - - - - - -
2295	// read photons and "map" them into the pixels
2296	//--------------------------------------------------
2297
2298	// initialize CPhoton
2299
2300	photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
2301
2302	// read the photons data
2303
2304	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2305
2306	// loop over the photons
2307
2308	while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
2309
2310	memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
2311
2312	fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
2313
2314	// increase number of photons
2315
2316	(*incph)++;
2317
2318	//
2319	// Pixelization
2320	//
2321
2322	cx = photon.get_x();
2323	cy = photon.get_y();
2324
2325	// get wavelength
2326
2327	wl = photon.get_wl();
2328
2329	// cout << "wl " << wl << " radius " << sqrt(cxcx + cycy) << "\n";
2330
2331	// check if photon has valid wavelength and is inside outermost camera radius
2332
2333	if( (wl > maxwl_nm) \|\| (wl < minwl_nm) \|\| (sqrt(cxcx + cycy) > radius ) ){
2334
2335	// cout << " lost first check\n";
2336
2337	// read next CPhoton
2338
2339	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2340
2341	// go to beginning of loop, the photon is lost
2342	continue;
2343
2344	}
2345
2346	ipixnum = -1;
2347
2348	for(i=0; i<cam->inumpixels; i++){
2349	if( bpoint_is_in_pix( cx, cy, i, cam) ){
2350	ipixnum = i;
2351	break;
2352	}
2353	}
2354
2355	if(ipixnum==-1){// the photon is in none of the pixels
2356
2357	// cout << " lost pixlization\n";
2358
2359	// read next CPhoton
2360
2361	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2362
2363	// go to beginning of loop, the photon is lost
2364	continue;
2365	}
2366
2367	//+++
2368	// QE simulation
2369	//---
2370
2371	// set pointer to the QE table of the relevant pixel
2372
2373	qept = (float **)QE[ipixnum];
2374
2375	// check if wl is inside table; outside the table, QE is assumed to be zero
2376
2377	if((wl < qept[0][0]) \|\| (wl > qept[0][pointsQE-1])){
2378
2379	// cout << " lost\n";
2380
2381	// read next Photon
2382
2383	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2384
2385	// go to beginning of loop
2386	continue;
2387
2388	}
2389
2390	// find data point in the QE table (-> k)
2391
2392	k = 1; // start at 1 because the condition was already tested for 0
2393	while (k < pointsQE-1 && qept[0][k] < wl){
2394	k++;
2395	}
2396
2397	// calculate the qe value between 0. and 1.
2398
2399	qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
2400
2401	// if random > quantum efficiency, reject it
2402
2403	if ( RandomNumber > qe ) {
2404
2405	// cout << " lost\n";
2406
2407	// read next Photon
2408
2409	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2410
2411	// go to beginning of loop
2412	continue;
2413
2414	}
2415
2416	//+++
2417	// The photon has produced a photo electron
2418	//---
2419
2420	// cout << " accepted\n";
2421
2422	// increment the number of photoelectrons in the relevant pixel
2423
2424	nphe[ipixnum] += 1.0;
2425
2426	t = photon.get_t() ;
2427
2428	// cout << " t " << t;
2429
2430	// find minimum and maximum arrival time
2431
2432	if(t < *tmin_ns){
2433	*tmin_ns = t; // memorize time
2434	}
2435	if(t > *tmax_ns){
2436	*tmax_ns = t; // memorize time
2437	}
2438
2439	// store the new photoelectron
2440
2441	if(*itotnphe >= iMAXNUMPHE){
2442	cout << "Error: Memory overflow. Event produces more than maximum\n";
2443	cout << " allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
2444	return(1);
2445	}
2446
2447	phe[*itotnphe].iarrtime_ns = (int)t;
2448	phe[*itotnphe].ipixnum = ipixnum;
2449
2450	*itotnphe += 1;
2451
2452	// read next Photon
2453
2454	fread( flag, SIZE_OF_FLAGS, 1, sp );
2455
2456	} // end while, i.e. found end of event
2457
2458	return(0);
2459
2460	}
2461
2462
2463	//------------------------------------------------------------
2464	// @name produce_nsbrates
2465	//
2466	// @desc read the starfield file, call produce_phes on it in,
2467	// @desc different wavebands, calculate the nsbrates
2468	//
2469	// @date Mon Feb 14 16:44:21 CET 2000
2470	// @function @code
2471	//------------------------------------------------------------
2472
2473	int produce_nsbrates( char *iname, // the starfield input file name
2474	struct camera *cam, // camera layout
2475	class Photoelectron phe[iMAXNUMPHE], // array for photoelectrons
2476	float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
2477	// the NSB rates in phe/ns for each pixel
2478	){
2479
2480	int i, j, k, ii; // counters
2481
2482	static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
2483	WAVEBANDBOUND2,
2484	WAVEBANDBOUND3,
2485	WAVEBANDBOUND4,
2486	WAVEBANDBOUND5,
2487	WAVEBANDBOUND6 };
2488
2489	FILE *infile; // the input file
2490	fpos_t fileposition; // marker on the input file
2491	static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
2492	static MCEventHeader evth; // the event header
2493	static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
2494	int itnphe; // total number of produced photoelectrons
2495	int itotnphe; // total number of produced photoelectrons after averaging
2496	int incph; // total number of cph read
2497	float tmin_ns; // minimum arrival time of all phes
2498	float tmax_ns; // maximum arrival time of all phes
2499	float integtime_ns; // integration time from the starfield generator
2500
2501	// open input file
2502
2503	log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
2504
2505	infile = fopen( iname, "r" );
2506	if ( infile == NULL )
2507	error( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
2508
2509	// get signature, and check it
2510
2511	if(check_reflector_file(infile)==FALSE){
2512	exit(1);
2513	}
2514
2515	// initialize flag
2516
2517	strcpy( cflag, " \0" );
2518
2519	// get flag
2520
2521	fread( cflag, SIZE_OF_FLAGS, 1, infile );
2522
2523	if ( ! feof(infile)){
2524
2525	// reading .rfl file
2526
2527	if(!isA( cflag, FLAG_START_OF_RUN )){
2528	error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
2529	}
2530	else { // found start of run
2531
2532	fread( cflag, SIZE_OF_FLAGS, 1, infile );
2533
2534	if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
2535
2536	// get MCEventHeader
2537
2538	fread( (char*)&evth, evth.mysize(), 1, infile );
2539
2540	integtime_ns = evth.get_energy();
2541
2542	// memorize where we are in the file
2543
2544	if (fgetpos( infile, &fileposition ) != 0){
2545	error( SIGNATURE, "Cannot position in file ...\n");
2546	}
2547
2548	// loop over the wavebands
2549
2550	for(i=0; i<iNUMWAVEBANDS; i++){
2551
2552	// initialize the rate array
2553
2554	for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
2555	rate_phepns[j][i] = 0.;
2556	}
2557
2558	itotnphe = 0;
2559
2560	// read the photons and produce the photoelectrons
2561	// - in order to average over the QE simulation, call the
2562	// production function iNUMNSBPRODCALLS times
2563
2564	for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
2565
2566	// position the file pointer to the beginning of the photons
2567
2568	fsetpos( infile, &fileposition);
2569
2570	// produce photoelectrons
2571
2572	k = produce_phes( infile,
2573	cam,
2574	wl_nm[i],
2575	wl_nm[i+1],
2576	phe, // this is a dummy here
2577	&itnphe,
2578	nphe, // we want this!
2579	&incph,
2580	&tmin_ns,
2581	&tmax_ns );
2582
2583	if( k != 0 ){ // non-zero returnvalue means error
2584	cout << "Exiting.\n";
2585	exit(1);
2586	}
2587
2588	for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
2589	rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
2590	}
2591
2592	itotnphe += itnphe;
2593
2594	} // end for(ii=0 ...
2595
2596	fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
2597	wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
2598
2599	} // end for(i=0 ...
2600
2601	}
2602	else{
2603	cout << "Starfield file contains no event.\nExiting.\n";
2604	exit(1);
2605	} // end if( isA ... event
2606	} // end if ( isA ... run
2607	}
2608	else{
2609	cout << "Starfield file contains no run.\nExiting.\n";
2610	exit(1);
2611	}
2612
2613	fclose( infile );
2614
2615	return(0);
2616
2617	}
2618
2619
2620	//------------------------------------------------------------
2621	// @name produce_nsb_phes
2622	//
2623	// @desc produce the photoelectrons from the nsbrates
2624	//
2625	// @date Thu Feb 17 17:10:40 CET 2000
2626	// @function @code
2627	//------------------------------------------------------------
2628
2629	int produce_nsb_phes( float *atmin_ns,
2630	float *atmax_ns,
2631	float theta_rad,
2632	struct camera *cam,
2633	float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
2634	float difnsbr_phepns[iMAXNUMPIX],
2635	float extinction[iNUMWAVEBANDS],
2636	float fnpx[iMAXNUMPIX],
2637	Photoelectron photo[iMAXNUMPHE],
2638	int *inphe,
2639	float base_mv[iMAXNUMPIX]){
2640
2641	float simtime_ns; // NSB simulation time
2642	int i, j, k, ii;
2643	float zenfactor; // correction factor calculated from the extinction
2644	int inumnsbphe; // number of photoelectrons caused by NSB
2645
2646	ii = *inphe; // avoid dereferencing
2647
2648	// check if the arrival times are set; if not generate them
2649
2650	if(*atmin_ns == NOTIME){
2651
2652	*atmin_ns = 0.;
2653	atmax_ns = simtime_ns = SLICESWIDTH_TIMESLICE;
2654
2655	}
2656	else{ // extend the event time window by the given offsets
2657
2658	atmin_ns = atmin_ns - SIMTIMEOFFSET_NS;
2659	atmax_ns = atmax_ns + SIMTIMEOFFSET_NS;
2660
2661	simtime_ns = atmax_ns - atmin_ns;
2662
2663	// make sure the simulated time is long enough for the FADC simulation
2664
2665	if(simtime_ns< SLICES*WIDTH_TIMESLICE){
2666	atmax_ns = atmin_ns + SLICES*WIDTH_TIMESLICE;
2667	simtime_ns = SLICES*WIDTH_TIMESLICE;
2668	}
2669
2670	}
2671
2672	// initialize baselines
2673
2674	for(i=0; i<cam->inumpixels; i++){
2675	base_mv[i] = 0.;
2676	}
2677
2678	// calculate baselines and generate phes
2679
2680	for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
2681
2682	// calculate the effect of the atmospheric extinction
2683
2684	zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
2685
2686	for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
2687
2688	inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
2689	* simtime_ns );
2690
2691	base_mv[j] += inumnsbphe;
2692
2693	// randomize
2694
2695	inumnsbphe = ignpoi( inumnsbphe );
2696
2697	// create the photoelectrons
2698
2699	for(k=0; k < inumnsbphe; k++){
2700
2701	if(ii >= iMAXNUMPHE){
2702	cout << "Error: Memory overflow. NSB simulation produces more than maximum\n";
2703	cout << " allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
2704	return(1);
2705	}
2706
2707	photo[ii].iarrtime_ns = (int)(RandomNumber * simtime_ns + *atmin_ns );
2708	photo[ii].ipixnum = j;
2709
2710	// cout << "Created phe " << photo[ii].iarrtime_ns << " "
2711	// << photo[ii].ipixnum << "\n";
2712
2713	ii++; // increment total number of photoelectons
2714
2715	fnpx[j] += 1.; // increment number of photoelectrons in this pixel
2716
2717	}
2718
2719	} // end for(j=0 ...
2720	} // end for(i=0 ...
2721
2722	// finish baseline calculation
2723
2724	for(i=0; i<cam->inumpixels; i++){
2725	base_mv[i] = RESPONSE_FWHM RESPONSE_AMPLITUDE / simtime_ns;
2726	}
2727
2728	*inphe = ii; // update the pointer
2729
2730	return(0);
2731	}
2732
2733
2734	// @endcode
2735
2736
2737	//=------------------------------------------------------------
2738	//!@subsection Log of this file.
2739
2740	//!@{
2741	//
2742	// $Log: not supported by cvs2svn $
2743	// Revision 1.6 2000/03/20 18:35:11 blanch
2744	// The trigger is already implemented but it does not save the trigger information in any file as it is implemented in timecam. In the next days there will be a version which also creates the files with the trigger information. It is going to be a mixing of the current camera and timecam programs.
2745	//
2746	// Revision 1.5 2000/02/18 17:40:35 petry
2747	// This version includes drastic changes compared to camera.cxx 1.4.
2748	// It is not yet finished and not immediately useful because the
2749	// trigger simulation is not yet re-implemented. I had to take it
2750	// out together with some other stuff in order to tidy the whole
2751	// program up. This is not meant as an insult to anyone. I needed
2752	// to do this in order to be able to work on it.
2753	//
2754	// This version has been put in the repository in order to be
2755	// able to share the further development with others.
2756	//
2757	// If you need something working, wait or take an earlier one.
2758	// See file README.
2759	//
2760	// Revision 1.4 2000/01/25 08:36:23 petry
2761	// The pixelization in previous versions was buggy.
2762	// This is the first version with a correct pixelization.
2763	//
2764	// Revision 1.3 2000/01/20 18:22:17 petry
2765	// Found little bug which makes camera crash if it finds a photon
2766	// of invalid wavelength. This bug is now fixed and the range
2767	// of valid wavelengths extended to 290 - 800 nm.
2768	// This is in preparation for the NSB simulation to come.
2769	// Dirk
2770	//
2771	// Revision 1.2 1999/11/19 08:40:42 harald
2772	// Now it is possible to compile the camera programm under osf1.
2773	//
2774	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
2775	// This the starting point for CVS controlled further developments of the
2776	// camera program. The program was originally written by Jose Carlos.
2777	// But here you can find a "rootified" version to the program. This means
2778	// that there is no hbook stuff in it now. Also the output of the
2779	// program changed to the MagicRawDataFormat.
2780	//
2781	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
2782	//
2783	// Revision 1.3 1999/10/22 15:01:28 petry
2784	// version sent to H.K. and N.M. on Fri Oct 22 1999
2785	//
2786	// Revision 1.2 1999/10/22 09:44:23 petry
2787	// first synthesized version which compiles and runs without crashing;
2788	//
2789	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
2790	// first synthesised version
2791	//
2792	// Revision 1.13 1999/03/15 14:59:05 gonzalez
2793	// camera-1_1
2794	//
2795	// Revision 1.12 1999/03/02 09:56:10 gonzalez
2796	// * empty log message *
2797	//
2798	//
2799	//!@}
2800
2801	//=EOF

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

Download in other formats: