Context Navigation

camera.cxx@ 373

Visit:

Last change on this file since 373 was 365, checked in by petry, 25 years ago

This version includes drastic changes compared to camera.cxx 1.4. It is not yet finished and not immediately useful because the trigger simulation is not yet re-implemented. I had to take it out together with some other stuff in order to tidy the whole program up. This is not meant as an insult to anyone. I needed to do this in order to be able to work on it. This version has been put in the repository in order to be able to share the further development with others. If you need something working, wait or take an earlier one. See file README.

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

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

Context Navigation

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

Download in other formats: