Context Navigation

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

Visit:

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

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

Download in other formats: