Context Navigation

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

Visit:

Last change on this file since 642 was 577, checked in by magicsol, 24 years ago
* empty log message *
File size: 86.0 KB

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

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

Download in other formats: