Context Navigation

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

Visit:

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

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

Download in other formats: