Context Navigation

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

Visit:

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

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

Download in other formats: