Context Navigation

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

Visit:

Last change on this file since 436 was 436, checked in by harald, 24 years ago
Added a lot of changes done by oscar.
File size: 83.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.12 $
24	// $Author: harald $
25	// $Date: 2000-09-22 17:40:18 $
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
745	if (Trigger_Loop && Write_McTrig){
746	for(char branchname[10],i=0;i<icontrigger;i++){
747	//
748	// build the name of the branch for the different trigger conditions
749	//
750	sprintf(help,"%i",i+1);
751	strcpy (branchname, "MHeaderTrig");
752	strcat (branchname, &help[0]);
753	strcat (branchname, ".");
754
755	HeaderTree.Branch(branchname,"MHeaderTrig",
756	&HeaderTrig[i], bsize, split);
757	}
758	}
759
760	// Fill branches
761
762	if(!Trigger_Loop && Write_McTrig){
763
764	HeaderTrig[0]->SetTopology((Short_t) Trigger_topology);
765	HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity);
766	for(i=0;i<TRIGGER_PIXELS;i++){
767	fpixelthres[i]=(Float_t)Trigger_threshold;
768	}
769	HeaderTrig[0]->SetThreshold( fpixelthres);
770
771	}
772
773
774	if(Trigger_Loop && Write_McTrig){
775	for (int iconcount=0,ithrescount=0;ithrescount<=Trigger_loop_uthres-Trigger_loop_lthres;ithrescount++){
776	for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
777	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
778	HeaderTrig[iconcount]->SetTopology((Short_t) itopocount+Trigger_loop_ltop);
779	HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
780	for(i=0;i<TRIGGER_PIXELS;i++){
781	fpixelthres[i]=(Float_t)(ithrescount+Trigger_loop_lthres);
782	}
783	HeaderTrig[iconcount]->SetThreshold( fpixelthres);
784	iconcount++;
785	}
786	}
787	}
788	}
789
790	// Fill the Header Tree with the current leaves of each branch
791
792	HeaderTree.Fill() ;
793
794	// create a Tree for the Event data stream
795	TTree EvtTree("EvtTree","Events of Run");
796
797	if (Write_McEvt){
798	EvtTree.Branch("MMcEvt","MMcEvt",
799	&McEvt, bsize, split);
800	}
801
802	if(!Trigger_Loop){
803	if (Write_RawEvt){
804	EvtTree.Branch("MRawEvt","MRawEvt",
805	&Evt[0], bsize, split);
806	}
807	if (Write_McTrig){
808	EvtTree.Branch("MMcTrig","MMcTrig",
809	&McTrig[0], bsize, split);
810	}
811	}
812	else{ // trigger lopp
813	if (Write_McTrig){
814	for(char branchname[10],i=0;i<icontrigger;i++){
815
816	sprintf(help,"%i",i+1);
817	strcpy (branchname, "MMcTrig");
818	strcat (branchname, & help[0]);
819	strcat (branchname, ".");
820	EvtTree.Branch(branchname,"MMcTrig",
821	&McTrig[i], bsize, split);
822	}
823	}
824	}
825
826	if (Trigger_Loop && Write_RawEvt){
827	for(char branchname[10],i=0;i<icontrigger;i++){
828
829	sprintf(help,"%i",i+1);
830	strcpy (branchname, "MRawEvt");
831	strcat (branchname, & help[0]);
832	strcat (branchname, ".");
833	EvtTree.Branch(branchname,"MRawEvt",
834	&Evt[i], bsize, split);
835	}
836	}
837
838	unsigned short ulli = 0 ;
839
840	TApplication theAppTrigger("App", &argc, argv);
841
842	if (gROOT->IsBatch()) {
843	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
844	// return 1;
845	}
846
847	if(FADC_Scan){
848	TApplication theAppFadc("App", &argc, argv);
849
850	if (gROOT->IsBatch()) {
851	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
852	// return 1;
853	}
854	}
855
856	// set plate scale (deg/cm) and trigger area (deg)
857
858	plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;
859
860	if ( ! get_trigger_radius( &degTriggerZone ) )
861	degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
862
863	if ( ! get_correction( &fCorrection ) )
864	fCorrection = 1.0;
865
866	// number of pixels for parameters
867
868	anaPixels = get_ana_pixels();
869	anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
870
871	// prepare the NSB simulation
872
873	if( simulateNSB ){
874
875	//
876	// Calculate the non-diffuse NSB photoelectron rates
877	//
878
879	k = produce_nsbrates( starfieldname,
880	&cam,
881	nsbrate_phepns );
882
883	if (k != 0){
884	cout << "Error when reading starfield... \nExiting.\n";
885	exit(1);
886	}
887
888	// calculate diffuse rate correcting for the pixel size
889
890	for(i=0; i<cam.inumpixels; i++){
891	diffnsb_phepns[i] = meanNSB *
892	cam.dpixsizefactor[i] * cam.dpixsizefactor[i];
893	}
894	}
895
896	//
897	// Read the reflector file with the Cherenkov data
898	//
899
900	// select input file
901
902	if ( Data_From_STDIN ) {
903	inputfile = stdin;
904	}
905	else{
906
907	log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
908	inputfile = fopen( inname, "r" );
909	if ( inputfile == NULL )
910	error( SIGNATURE, "Cannot open input file: %s\n", inname );
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	if ( fmod ( nshow, 1000. ) == 0. )
987	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
988
989	// get MCEventHeader
990
991	fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
992
993	// calculate core distance and impact parameter
994
995	coreD = mcevth.get_core(&coreX, &coreY);
996
997	// calculate impact parameter (shortest distance betwee the original
998	// trajectory of the primary (assumed shower-axis) and the
999	// direction where the telescope points to
1000	//
1001	// we use the following equation, given that the shower core position
1002	// is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
1003	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
1004	// and (l2,m2,n2)
1005	//
1006	// \| \|
1007	// \| x1-x2 y1-y2 z1-z2 \|
1008	// \| \|
1009	// + \| l1 m1 n1 \|
1010	// - \| \|
1011	// \| l2 m2 n2 \|
1012	// \| \|
1013	// dist = ------------------------------------ ( > 0 )
1014	// [ \|l1 m1\|2 \|m1 n1\|2 \|n1 l1\|2 ]1/2
1015	// [ \| \| + \| \| + \| \| ]
1016	// [ \|l2 m2\| \|m2 n2\| \|n2 l2\| ]
1017	//
1018	// playing a little bit, we get this reduced for in our case:
1019	//
1020	//
1021	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
1022	// [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
1023	// 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
1024
1025	// read the direction of the incoming shower
1026
1027	thetashw = mcevth.get_theta();
1028	phishw = mcevth.get_phi();
1029
1030	// calculate vector for shower
1031
1032	l1 = sin(thetashw)*cos(phishw);
1033	m1 = sin(thetashw)*sin(phishw);
1034	n1 = cos(thetashw);
1035
1036	// read the deviation of the telescope with respect to the shower
1037
1038	mcevth.get_deviations ( &thetaCT, &phiCT );
1039
1040	if ( (thetaCT == 0.) && (phiCT == 0.) ) {
1041
1042	// CT was looking to the source (both lines are parallel)
1043	// therefore, we calculate the impact parameter as the distance
1044	// between the CT axis and the core position
1045
1046	impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1047
1048	} else {
1049
1050	// the shower comes off-axis
1051
1052	// obtain with this the final direction of the CT
1053
1054	thetaCT += thetashw;
1055	phiCT += phishw;
1056
1057	// calculate vector for telescope
1058
1059	l2 = sin(thetaCT)*cos(phiCT);
1060	m2 = sin(thetaCT)*sin(phiCT);
1061	n2 = cos(thetaCT);
1062
1063	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1064	den = (SQR(l1m2 - l2m1) +
1065	SQR(m1n2 - m2n1) +
1066	SQR(n1l2 - n2l1));
1067	den = sqrt(den);
1068
1069	impactD = fabs(num)/den;
1070
1071	// fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ",
1072	// thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2,
1073	// num, den);
1074
1075	}
1076
1077	// energy cut
1078
1079	if ( Select_Energy ) {
1080	if (( mcevth.get_energy() < Select_Energy_le ) \|\|
1081	( mcevth.get_energy() > Select_Energy_ue )) {
1082	log(SIGNATURE, "select_energy: shower rejected.\n");
1083	continue;
1084	}
1085	}
1086
1087	// Read first and last time and put inumphe to 0
1088
1089	mcevth.get_times(&arrtmin_ns,&arrtmax_ns);
1090	inumphe=0;
1091
1092	// NSB simulation
1093
1094	if(simulateNSB){
1095
1096	k = produce_nsb_phes( &arrtmin_ns, // will be changed by the function!
1097	&arrtmax_ns, // will be changed by the function!
1098	thetaCT,
1099	&cam,
1100	nsbrate_phepns,
1101	diffnsb_phepns,
1102	ext,
1103	fnpix, // will be changed by the function!
1104	&Trigger, // will be changed by the function!
1105	&fadc, // will be changed by the function!
1106	&inumphe, // important for later: the size of photoe[]
1107	baseline_mv // will be generated by the function
1108	);
1109
1110	if( k != 0 ){ // non-zero returnvalue means error
1111	cout << "Exiting.\n";
1112	exit(1);
1113	}
1114
1115	}// end if(simulateNSB) ...
1116
1117	// read the photons and produce the photoelectrons
1118
1119	k = produce_phes( inputfile,
1120	&cam,
1121	WAVEBANDBOUND1,
1122	WAVEBANDBOUND6,
1123	&Trigger, // will be changed by the function!
1124	&fadc, // will be changed by the function!
1125	&inumphe, // important for later: the size of photoe[]
1126	fnpix, // will be changed by the function!
1127	&ncph, // will be changed by the function!
1128	&arrtmin_ns, // will be changed by the function!
1129	&arrtmax_ns // will be changed by the function!
1130	);
1131
1132	if( k != 0 ){ // non-zero returnvalue means error
1133	cout << "Exiting.\n";
1134	exit(1);
1135	}
1136
1137	if ( fmod ( nshow, 1000. ) == 0. )
1138	log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1139	ncph, ntcph);
1140
1141	ntcph += ncph;
1142
1143	// skip it ?
1144
1145	for ( i=0; i<nSkip; ++i ) {
1146	if (Skip[i] == (nshow+ntshow)) {
1147	i = -1;
1148	break;
1149	}
1150	}
1151
1152	// if after the previous loop, the exit value of i is -1
1153	// then the shower number is in the list of showers to be
1154	// skipped
1155
1156	if (i == -1) {
1157	log(SIGNATURE, "\t\tskipped!\n");
1158	continue;
1159	}
1160
1161	if ( fmod ( nshow, 1000. ) == 0. )
1162	log(SIGNATURE, "Total number of phe: %d \n", inumphe ) ;
1163
1164
1165	//++++++++++++++++++++++++++++++++++++++++++++++++++
1166	// at this point we have a camera full of
1167	// ph.e.s
1168	// we should first apply the trigger condition,
1169	// and if there's trigger, then clean the image,
1170	// calculate the islands statistics and the
1171	// other parameters of the image (Hillas' parameters
1172	// and so on).
1173	//--------------------------------------------------
1174
1175	// TRIGGER HERE
1176
1177	//
1178	// now the noise of the electronic
1179	// (preamps, optical transmission,..) is introduced.
1180	// This is done inside the class MTrigger by the method ElecNoise.
1181	//
1182	Trigger.ElecNoise() ;
1183
1184	fadc.ElecNoise() ;
1185
1186	// We study several trigger conditons
1187	if(Trigger_Loop){
1188	// Loop over trigger threshold
1189	for (int iconcount=0,ithrescount=Trigger_loop_lthres;ithrescount<=Trigger_loop_uthres;ithrescount++) {
1190	for (i=0;i<TRIGGER_PIXELS;i++)
1191	fpixelthres[i]=(float) ithrescount;
1192
1193	Trigger.SetThreshold(fpixelthres);
1194	Trigger.Diskriminate();
1195
1196	//
1197	// look if in all the signals in the trigger signal branch
1198	// is a possible Trigger. Therefore we habe to diskriminate all
1199	// the simulated analog signals (Method Diskriminate in class
1200	// MTrigger). We look simultanously for the moments at which
1201	// there are more than TRIGGER_MULTI pixels above the
1202	// CHANNEL_THRESHOLD.
1203	//
1204
1205	// Set trigger flags to zero
1206	Lev1=Lev2=0;
1207	btrigger=0;
1208
1209	// loop over multiplicity of trigger configuration
1210	for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1211	Trigger.SetMultiplicity(imulticount);
1212	Trigger.ClearZero();
1213
1214	Lev0=(Short_t) Trigger.ZeroLevel();
1215	if (Lev0>0 \|\| Write_All_Images \|\| btrigger){
1216
1217	// loop over topologies
1218	for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1219	Lev1=Lev2=0;
1220
1221	if(itopocount==0 && imulticount>7) continue;
1222	if(itopocount==2 && imulticount<3) continue;
1223	Trigger.SetTopology(itopocount);
1224	Trigger.ClearFirst();
1225
1226	//
1227	// Start the First Level Trigger simulation
1228	//
1229	Lev1=Trigger.FirstLevel();
1230	if (Write_McTrig)
1231	McTrig[iconcount]->SetFirstLevel (Lev1);
1232	if(Lev1>0) {
1233	btrigger= 1;
1234	ntriggerloop[ithrescount-Trigger_loop_lthres][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++;
1235	}
1236
1237	if(Lev1==0 && (Write_All_Images \|\| btrigger)){
1238	btrigger= 1;
1239	Lev1=1;
1240	}
1241	numPix=0;
1242	for (Int_t ii=0;ii<Lev1;ii++){
1243	if (Write_McTrig)
1244	McTrig[iconcount]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
1245	if (Write_McTrig)
1246	McTrig[iconcount]->SetPixel(Trigger.GetFirstLevelPixel(ii),ii+1);
1247
1248	//
1249	// fill inside the class fadc the member output
1250	//
1251
1252	fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
1253
1254	if( Write_RawEvt ){
1255	//
1256	// Fill the header of this event
1257	//
1258
1259	Evt[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),(Float_t) (nshow10000+iconcount100+ii),0 );
1260	// fill pixel information
1261	for(i=0;i<ct_NPixels;i++){
1262	for (j=0;j<SLICES_MFADC;j++){
1263	fadcValues[j]=fadc.GetFadcSignal(i,j);
1264	}
1265	Evt[iconcount]->FillPixel(i+1000*ii,&numPix,fadcValues);
1266	}
1267	}
1268	}
1269	//
1270	// Increase counter of analised trigger conditions
1271	//
1272	iconcount++;
1273	}
1274	}
1275	else{
1276	break;
1277	}
1278	}
1279	if (!btrigger) break;
1280	}
1281	if (btrigger){
1282
1283	//
1284	// fill the MMcEvt with all information
1285	//
1286
1287	if (Write_McEvt) {
1288	McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1289	mcevth.get_energy(),
1290	mcevth.get_theta(),
1291	mcevth.get_phi(),
1292	mcevth.get_core(),
1293	mcevth.get_coreX(),
1294	mcevth.get_coreY(),
1295	impactD,
1296	ulli, ulli,
1297	(UShort_t) ncph,
1298	ulli,
1299	(UShort_t) ncph) ;
1300	}
1301	// Fill the Tree with the current leaves of each branch
1302	i=EvtTree.Fill() ;
1303
1304	// Clear the branches
1305	if(Write_McTrig){
1306	for(i=0;i<icontrigger;i++){
1307	McTrig[i]->Clear() ;
1308	}
1309	}
1310	if( Write_RawEvt ){
1311	for(i=0;i<icontrigger;i++){
1312	Evt[i]->Clear() ;
1313	}
1314	}
1315	if (Write_McEvt)
1316	McEvt->Clear() ;
1317	}
1318	}
1319	// We study a single trigger condition
1320	else {
1321
1322	// Setting trigger conditions
1323	Trigger.SetMultiplicity(Trigger_multiplicity);
1324	Trigger.SetTopology(Trigger_topology);
1325	for (i=0;i<TRIGGER_PIXELS;i++)
1326	fpixelthres[i]=Trigger_threshold;
1327	Trigger.SetThreshold(fpixelthres);
1328
1329	Trigger.Diskriminate() ;
1330
1331	//
1332	// look if in all the signals in the trigger signal branch
1333	// is a possible Trigger. Therefore we habe to diskriminate all
1334	// the simulated analog signals (Method Diskriminate in class
1335	// MTrigger). We look simultanously for the moments at which
1336	// there are more than TRIGGER_MULTI pixels above the
1337	// CHANNEL_THRESHOLD.
1338	//
1339
1340	Lev0 = (Short_t) Trigger.ZeroLevel() ;
1341
1342	Lev1 = Lev2 = 0 ;
1343
1344	//
1345	// Start the First Level Trigger simulation
1346	//
1347
1348	if ( Lev0 > 0 \|\| Write_All_Images) {
1349	Lev1= Trigger.FirstLevel();
1350	if (Write_McTrig)
1351	McTrig[0]->SetFirstLevel (Lev1);
1352	}
1353	if (Lev1>0){
1354	++ntrigger;
1355	}
1356	if (Lev1==0 && Write_All_Images){
1357	Lev1=1;
1358	}
1359
1360	numPix=0;
1361	for(Int_t ii=0;ii<Lev1;ii++){
1362	// Loop over different level one triggers
1363
1364	//
1365	// fill inside class fadc the member output
1366	//
1367	fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
1368
1369	if (Write_McTrig)
1370	McTrig[0]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
1371
1372	if (Write_McTrig)
1373	McTrig[0]->SetPixel(Trigger.GetFirstLevelPixel(ii),ii+1);
1374
1375	// Fill Evt information
1376
1377	if (Write_RawEvt){
1378
1379	//
1380	// Fill the header of this event
1381	//
1382
1383	Evt[0]->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ;
1384
1385	// fill pixel information
1386
1387	for(i=0;i<ct_NPixels;i++){
1388	for (j=0;j<SLICES_MFADC;j++){
1389	fadcValues[j]=fadc.GetFadcSignal(i,j);
1390	}
1391	Evt[0]->FillPixel(i,&numPix,fadcValues);
1392	}
1393	}
1394	//
1395	// fill the MMcEvt with all information
1396	//
1397
1398	if (Write_McEvt){
1399	McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1400	mcevth.get_energy(),
1401	mcevth.get_theta(),
1402	mcevth.get_phi(),
1403	mcevth.get_core(),
1404	mcevth.get_coreX(),
1405	mcevth.get_coreY(),
1406	impactD,
1407	ulli, ulli,
1408	(UShort_t) ncph,
1409	ulli,
1410	(UShort_t) ncph) ;
1411	}
1412
1413	// We don not count photons out of the camera.
1414
1415	//
1416	// write it out to the file outfile
1417	//
1418
1419	//EvtTree.Fill() ;
1420	// huschel
1421
1422
1423	//
1424	// if a first level trigger occurred, then
1425	// 1. do some other stuff (not implemented)
1426	// 2. start the gui tool
1427
1428	if(FADC_Scan){
1429	if ( Lev0 > 0 ) {
1430	fadc.ShowSignal( McEvt, (Float_t) 60. ) ;
1431	}
1432	}
1433
1434	if(Trigger_Scan){
1435	if ( Lev0 > 0 ) {
1436	Trigger.ShowSignal(McEvt) ;
1437	}
1438	}
1439
1440	// clear all
1441	if (Write_RawEvt) Evt[0]->Clear() ;
1442	if (Write_McEvt) McEvt->Clear() ;
1443	}
1444	if (Write_McTrig) McTrig[0]->Clear() ;
1445	}
1446
1447	#ifdef __DEBUG__
1448	printf("\n");
1449
1450	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1451
1452	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1453
1454	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1455
1456	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1457
1458	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1459	(int)pixels[ici][icj][PIXNUM],
1460	pixels[ici][icj][PIXX],
1461	pixels[ici][icj][PIXY],
1462	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1463
1464	}
1465	}
1466	}
1467	}
1468
1469	for (i=0; i<ct_NPixels; ++i) {
1470	printf("%d (%d): ", i, npixneig[i]);
1471	for (j=0; j<npixneig[i]; ++i)
1472	printf(" %d", pixneig[i][j]);
1473	printf("\n");
1474	}
1475
1476	#endif // __DEBUG__
1477
1478
1479	// look for the next event
1480
1481	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1482
1483	} // end while there is a next event
1484
1485	if( !isA( flag, FLAG_END_OF_RUN )){
1486	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1487	}
1488	else { // found end of run
1489	ntshow += nshow;
1490	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1491
1492	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1493
1494	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1495	log(SIGNATURE, "End of file . . .\n");
1496	still_in_loop = FALSE;
1497
1498	if ((! Data_From_STDIN) && ( !feof(inputfile) )){
1499
1500	// we have concatenated input files.
1501	// get signature of the next part and check it.
1502
1503	if(check_reflector_file( inputfile )==FALSE){
1504	exit(1);
1505	}
1506
1507	}
1508
1509	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1510
1511	} // end if found end of file
1512	} // end if found end of run
1513
1514	} // end if else found start of run
1515	} // end big while loop
1516
1517	//++
1518	// put the Event to the root file
1519	//--
1520
1521	HeaderTree.Write() ;
1522	EvtTree.Write() ;
1523	outfile_temp.Write() ;
1524	outfile_temp.Close() ;
1525
1526	//
1527	// initalize the ROOT file
1528	//
1529
1530	TFile outfile ( rootname , "RECREATE" );
1531	HeaderTree.Write() ;
1532	EvtTree.Write() ;
1533	outfile.Write() ;
1534	outfile.Close() ;
1535
1536	// close input file
1537
1538	log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
1539	ntshow, ntcph );
1540	datafile<<ntshow<<" event(s), with a total of "<<ntcph<<" C.photons"<<endl;
1541	if (Trigger_Loop){
1542	log( SIGNATURE, "Trigger Mode. \n");
1543	log( SIGNATURE, "Fraction of triggers: \n");
1544	datafile<<"Fraction of triggers: "<<endl;
1545	for (ithrescount=Trigger_loop_lthres;ithrescount<=Trigger_loop_uthres;ithrescount++){
1546	for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1547	for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1548	log( SIGNATURE, "Thres %d, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
1549	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);
1550	datafile<<"Thres "<<ithrescount<<", Multi "<<imulticount<<", Topo"<<itopocount<<": ";
1551	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;
1552	}
1553	}
1554	}
1555	}
1556	else{
1557	log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
1558	((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
1559	datafile<<"Fraction of triggers: "<<((float)ntrigger) / ((float)ntshow) * 100.0<<" ("<<ntrigger<<" out of "<<ntshow<<" )"<<endl;
1560	}
1561
1562	// close files
1563
1564	log( SIGNATURE, "Closing files\n" );
1565
1566	if( ! Data_From_STDIN ){
1567	fclose( inputfile );
1568	}
1569	datafile.close();
1570
1571	// program finished
1572
1573	log( SIGNATURE, "Done.\n");
1574
1575	return( 0 );
1576	}
1577	//!@}
1578
1579	// @T \newpage
1580
1581	//!@subsection Functions definition.
1582
1583	//!-----------------------------------------------------------
1584	// @name present
1585	//
1586	// @desc Make some presentation
1587	//
1588	// @date Sat Jun 27 05:58:56 MET DST 1998
1589	//------------------------------------------------------------
1590	// @function
1591
1592	//!@{
1593	void
1594	present(void)
1595	{
1596	cout << "##################################################\n"
1597	<< SIGNATURE << '\n' << '\n'
1598	<< "Processor of the reflector output\n"
1599	<< "J C Gonzalez, Jun 1998\n"
1600	<< "##################################################\n\n"
1601	<< flush ;
1602	}
1603	//!@}
1604
1605
1606	//!-----------------------------------------------------------
1607	// @name usage
1608	//
1609	// @desc show help
1610	//
1611	// @date Tue Dec 15 16:23:30 MET 1998
1612	//------------------------------------------------------------
1613	// @function
1614
1615	//!@{
1616	void
1617	usage(void)
1618	{
1619	present();
1620	cout << "\nusage ::\n\n"
1621	<< "\t camera "
1622	<< " [ -@ paramfile ] "
1623	<< " [ -h ] "
1624	<< "\n\n or \n\n"
1625	<< "\t camera < paramfile"
1626	<< "\n\n";
1627	exit(0);
1628	}
1629	//!@}
1630
1631
1632	//!-----------------------------------------------------------
1633	// @name log
1634	//
1635	// @desc function to send log information
1636	//
1637	// @var funct Name of the caller function
1638	// @var fmt Format to be used (message)
1639	// @var ... Other information to be shown
1640	//
1641	// @date Sat Jun 27 05:58:56 MET DST 1998
1642	//------------------------------------------------------------
1643	// @function
1644
1645	//!@{
1646	void
1647	log(const char funct, char fmt, ...)
1648	{
1649	va_list args;
1650
1651	// Display the name of the function that called error
1652	printf("[%s]: ", funct);
1653
1654	// Display the remainder of the message
1655	va_start(args, fmt);
1656	vprintf(fmt, args);
1657	va_end(args);
1658	}
1659	//!@}
1660
1661
1662	//!-----------------------------------------------------------
1663	// @name error
1664	//
1665	// @desc function to send an error message, and abort the program
1666	//
1667	// @var funct Name of the caller function
1668	// @var fmt Format to be used (message)
1669	// @var ... Other information to be shown
1670	//
1671	// @date Sat Jun 27 05:58:56 MET DST 1998
1672	//------------------------------------------------------------
1673	// @function
1674
1675	//!@{
1676	void
1677	error(const char funct, char fmt, ...)
1678	{
1679	va_list args;
1680
1681	// Display the name of the function that called error
1682	fprintf(stdout, "ERROR in %s: ", funct);
1683
1684	// Display the remainder of the message
1685	va_start(args, fmt);
1686	vfprintf(stdout, fmt, args);
1687	va_end(args);
1688
1689	perror(funct);
1690
1691	exit(1);
1692	}
1693	//!@}
1694
1695
1696	//!-----------------------------------------------------------
1697	// @name isA
1698	//
1699	// @desc returns TRUE(FALSE), if the flag is(is not) the given
1700	//
1701	// @var s1 String to be searched
1702	// @var flag Flag to compare with string s1
1703	// @return TRUE: both strings match; FALSE: oth.
1704	//
1705	// @date Wed Jul 8 15:25:39 MET DST 1998
1706	//------------------------------------------------------------
1707	// @function
1708
1709	//!@{
1710	int
1711	isA( char * s1, const char * flag ) {
1712	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
1713	}
1714	//!@}
1715
1716
1717	//!-----------------------------------------------------------
1718	// @name read_ct_file
1719	//
1720	// @desc read CT definition file
1721	//
1722	// @date Sat Jun 27 05:58:56 MET DST 1998
1723	//------------------------------------------------------------
1724	// @function
1725
1726	//!@{
1727	void
1728	read_ct_file(void)
1729	{
1730	char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
1731	char token[ITEM_MAX_LENGTH]; //@< a single token
1732	int i, j; //@< dummy counters
1733
1734	log( "read_ct_file", "start.\n" );
1735
1736	ifstream ctin ( ct_filename );
1737
1738	if ( ctin.bad() )
1739	error( "read_ct_file",
1740	"Cannot open CT def. file: %s\n", ct_filename );
1741
1742	// loop till the "end" directive is reached
1743
1744	while (!ctin.eof()) {
1745
1746	// get line from stdin
1747
1748	ctin.getline(line, LINE_MAX_LENGTH);
1749
1750	// look for each item at the beginning of the line
1751
1752	for (i=0; i<=define_mirrors; i++)
1753	if (strstr(line, CT_ITEM_NAMES[i]) == line)
1754	break;
1755
1756	// if it is not a valid line, just ignore it
1757
1758	if (i == define_mirrors+1)
1759	continue;
1760
1761	// case block for each directive
1762
1763	switch ( i ) {
1764
1765	case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
1766
1767	// get focal distance
1768
1769	sscanf(line, "%s %d", token, &ct_Type);
1770
1771	log( "read_ct_file", "<Type of Telescope>: %s\n",
1772	((ct_Type==0) ? "CT1" : "MAGIC") );
1773
1774	break;
1775
1776	case focal_distance: // <focal distance> [cm]
1777
1778	// get focal distance
1779
1780	sscanf(line, "%s %f", token, &ct_Focal_mean);
1781
1782	log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
1783
1784	break;
1785
1786	case focal_std: // s(focal distance) [cm]
1787
1788	// get focal distance
1789
1790	sscanf(line, "%s %f", token, &ct_Focal_std);
1791
1792	log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
1793
1794	break;
1795
1796	case point_spread: // <point spread> [cm]
1797
1798	// get point spread
1799
1800	sscanf(line, "%s %f", token, &ct_PSpread_mean);
1801
1802	log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
1803
1804	break;
1805
1806	case point_std: // s(point spread) [cm]
1807
1808	// get point spread
1809
1810	sscanf(line, "%s %f", token, &ct_PSpread_std);
1811
1812	log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
1813
1814	break;
1815
1816	case adjustment_dev: // s(adjustment_dev) [cm]
1817
1818	// get point spread
1819
1820	sscanf(line, "%s %f", token, &ct_Adjustment_std);
1821
1822	log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
1823
1824	break;
1825
1826	case black_spot: // radius of the black spot in the center [cm]
1827
1828	// get black spot radius
1829
1830	sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
1831
1832	log( "read_ct_file", "Radius of the black spots: %f cm\n",
1833	ct_BlackSpot_rad);
1834
1835	break;
1836
1837	case r_mirror: // radius of the mirrors [cm]
1838
1839	// get radius of mirror
1840
1841	sscanf(line, "%s %f", token, &ct_RMirror);
1842
1843	log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
1844
1845	break;
1846
1847	case n_mirrors: // number of mirrors
1848
1849	// get the name of the output_file from the line
1850
1851	sscanf(line, "%s %d", token, &ct_NMirrors);
1852
1853	log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
1854
1855	break;
1856
1857	case camera_width: // camera width [cm]
1858
1859	// get the name of the ct_file from the line
1860
1861	sscanf(line, "%s %f", token, &ct_CameraWidth);
1862
1863	log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
1864
1865	break;
1866
1867	case n_pixels: // number of pixels
1868
1869	// get the name of the output_file from the line
1870
1871	sscanf(line, "%s %d", token, &ct_NPixels);
1872
1873	log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
1874
1875	break;
1876
1877	case n_centralpixels: // number of central pixels
1878
1879	// get the name of the output_file from the line
1880
1881	sscanf(line, "%s %d", token, &ct_NCentralPixels);
1882
1883	log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
1884
1885	break;
1886
1887	case n_gappixels: // number of gap pixels
1888
1889	// get the name of the output_file from the line
1890
1891	sscanf(line, "%s %d", token, &ct_NGapPixels);
1892
1893	log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
1894
1895	break;
1896
1897	case pixel_width: // pixel width [cm]
1898
1899	// get the name of the ct_file from the line
1900
1901	sscanf(line, "%s %f", token, &ct_PixelWidth);
1902
1903	ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
1904	ct_PixelWidth_corner_2_corner_half =
1905	ct_PixelWidth_corner_2_corner * 0.50;
1906	ct_Apot = ct_PixelWidth / 2;
1907	ct_2Apot = ct_Apot * 2.0;
1908
1909	log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
1910
1911	break;
1912
1913	case define_mirrors: // read table with the parameters of the mirrors
1914
1915	log( "read_ct_file", "Table of mirrors data:\n" );
1916
1917	// check whether the number of mirrors was already set
1918
1919	if ( ct_NMirrors == 0 )
1920	error( "read_ct_file", "NMirrors was not set.\n" );
1921
1922	// allocate memory for paths list
1923
1924	log( "read_ct_file", "Allocating memory for ct_data\n" );
1925
1926	ct_data = new float*[ct_NMirrors];
1927
1928	for (i=0; i<ct_NMirrors; i++)
1929	ct_data[i] = new float[CT_NDATA];
1930
1931	// read data
1932
1933	log( "read_ct_file", "Reading mirrors data...\n" );
1934
1935	for (i=0; i<ct_NMirrors; i++)
1936	for (j=0; j<CT_NDATA; j++)
1937	ctin >> ct_data[i][j];
1938
1939	break;
1940
1941	} // switch ( i )
1942
1943	} // end while
1944
1945	// end
1946
1947	log( "read_ct_file", "done.\n" );
1948
1949	return;
1950	}
1951	//!@}
1952
1953
1954	//!-----------------------------------------------------------
1955	// @name read_pixels
1956	//
1957	// @desc read pixels data
1958	//
1959	// @date Fri Mar 12 16:33:34 MET 1999
1960	//------------------------------------------------------------
1961	// @function
1962
1963	//!@{
1964	void
1965	read_pixels(struct camera *pcam)
1966	{
1967	ifstream qefile;
1968	char line[LINE_MAX_LENGTH];
1969	int n, i, j, icount;
1970	float qe;
1971
1972	//------------------------------------------------------------
1973	// first, pixels' coordinates
1974
1975	pcam->inumpixels = ct_NPixels;
1976	pcam->inumcentralpixels = ct_NCentralPixels;
1977	pcam->inumgappixels = ct_NGapPixels;
1978	pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
1979	pcam->dpixdiameter_cm = ct_PixelWidth;
1980
1981	// initialize pixel numbers
1982
1983	pixary = new float* [2*ct_NCentralPixels];
1984	for ( i=0; i<2*ct_NCentralPixels; ++i )
1985	pixary[i] = new float[2];
1986
1987	pixneig = new int* [ct_NCentralPixels];
1988	for ( i=0; i<ct_NCentralPixels; ++i ) {
1989	pixneig[i] = new int[6];
1990	for ( j=0; j<6; ++j )
1991	pixneig[i][j] = -1;
1992	}
1993
1994	npixneig = new int[ct_NCentralPixels];
1995	for ( i=0; i<ct_NCentralPixels; ++i )
1996	npixneig[i] = 0;
1997
1998	// generate all coordinates
1999
2000	igen_pixel_coordinates(pcam);
2001
2002
2003	// calculate tables of neighbours
2004
2005	#ifdef __DEBUG__
2006	for ( n=0 ; n<ct_NPixels ; ++n ) {
2007	cout << "Para el pixel " << n << ": ";
2008	for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2009	if ( pixels_are_neig(n,i) == TRUE ) {
2010	pixneig[n][npixneig[n]] = i;
2011	pixneig[i][npixneig[i]] = n;
2012	cout << i << ' ';
2013	++npixneig[n];
2014	++npixneig[i];
2015	}
2016	}
2017	cout << endl << flush;
2018	}
2019	#else // ! __DEBUG__
2020	for ( n=0 ; n<ct_NCentralPixels ; ++n )
2021	for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2022	if ( pixels_are_neig(n,i) == TRUE ) {
2023	pixneig[n][npixneig[n]] = i;
2024	pixneig[i][npixneig[i]] = n;
2025	++npixneig[n];
2026	++npixneig[i];
2027	}
2028	#endif // ! __DEBUG__
2029
2030	#ifdef __DEBUG__
2031	for ( n=0 ; n<ct_NPixels ; ++n ) {
2032	cout << n << ':';
2033	for ( j=0; j<npixneig[n]; ++j)
2034	cout << ' ' << pixneig[n][j];
2035	cout << endl << flush;
2036	}
2037	#endif // __DEBUG__
2038
2039	//------------------------------------------------------------
2040	// second, pixels' QE
2041
2042	// try to open the file
2043
2044	log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
2045
2046	qefile.open( QE_FILE );
2047
2048	// if it is wrong or does not exist, exit
2049
2050	if ( qefile.bad() )
2051	error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2052
2053	// read file
2054
2055	log("read_pixels", "Reading data . . .\n");
2056
2057	i=-1;
2058	icount = 0;
2059
2060	while ( ! qefile.eof() ) {
2061
2062	// get line from the file
2063
2064	qefile.getline(line, LINE_MAX_LENGTH);
2065
2066	// skip if comment
2067
2068	if ( *line == '#' )
2069	continue;
2070
2071	// if it is the first valid value, it is the number of QE data points
2072
2073	if ( i < 0 ) {
2074
2075	// get the number of datapoints
2076
2077	sscanf(line, "%d", &pointsQE);
2078
2079	// allocate memory for the table of QEs
2080
2081	QE = new float ** [ct_NPixels];
2082
2083	for ( i=0; i<ct_NPixels; ++i ) {
2084	QE[i] = new float * [2];
2085	QE[i][0] = new float[pointsQE];
2086	QE[i][1] = new float[pointsQE];
2087	}
2088
2089	QElambda = new float [pointsQE];
2090
2091	for ( i=0; i<pointsQE; ++i ) {
2092	qefile.getline(line, LINE_MAX_LENGTH);
2093	sscanf(line, "%f", &QElambda[i]);
2094	}
2095
2096	i=0;
2097
2098	continue;
2099	}
2100
2101	// get the values (num-pixel, num-datapoint, QE-value)
2102
2103	if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2104	break;
2105
2106	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2107	((j-1) < pointsQE) && ((j-1) > -1) ) {
2108	QE[i-1][0][j-1] = QElambda[j-1];
2109	QE[i-1][1][j-1] = qe;
2110	}
2111
2112	if ( i > ct_NPixels) break;
2113
2114	icount++;
2115
2116	}
2117
2118	if(icount/pointsQE < ct_NPixels){
2119	error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2120	icount/pointsQE, ct_NPixels );
2121	}
2122
2123	// close file
2124
2125	qefile.close();
2126
2127	// test QE
2128
2129	for(icount=0; icount< ct_NPixels; icount++){
2130	for(i=0; i<pointsQE; i++){
2131	if( QE[icount][0][i] < 100. \|\| QE[icount][0][i] > 1000. \|\|
2132	QE[icount][1][i] < 0. \|\| QE[icount][1][i] > 100.){
2133	error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2134	icount, i, QE[icount][0][i], QE[icount][1][i] );
2135	}
2136	}
2137	}
2138
2139	// end
2140
2141	log("read_pixels", "Done.\n");
2142
2143	}
2144	//!@}
2145
2146
2147	//!-----------------------------------------------------------
2148	// @name pixels_are_neig
2149	//
2150	// @desc check whether two pixels are neighbours
2151	//
2152	// @var pix1 Number of the first pixel
2153	// @var pix2 Number of the second pixel
2154	// @return TRUE: both pixels are neighbours; FALSE: oth.
2155	//
2156	// @date Wed Sep 9 17:58:37 MET DST 1998
2157	//------------------------------------------------------------
2158	// @function
2159
2160	//!@{
2161	int
2162	pixels_are_neig(int pix1, int pix2)
2163	{
2164	if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2165	SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2166	> ct_PixelWidth_corner_2_corner )
2167	return ( FALSE );
2168	else
2169	return ( TRUE );
2170	}
2171	//!@}
2172
2173	//!-----------------------------------------------------------
2174	// @name igen_pixel_coordinates
2175	//
2176	// @desc generate the pixel center coordinates
2177	//
2178	// @var *pcam structure camera containing all the
2179	// camera information
2180	// @return total number of pixels
2181	//
2182	// DP
2183	//
2184	// @date Thu Oct 14 10:41:03 CEST 1999
2185	//------------------------------------------------------------
2186	// @function
2187
2188	//!@{
2189	/****** igen_pixel_coordinates() *******************************/
2190
2191	int igen_pixel_coordinates(struct camera *pcam) {
2192	/* generate pixel coordinates, return value is number of pixels */
2193
2194	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2195	float fsegment_fract;
2196	double dtsize;
2197	double dhsize;
2198	double dpsize;
2199	double dxfirst_pix;
2200	double dyfirst_pix;
2201	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
2202	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
2203
2204
2205	double dstartx, dstarty; /* for the gap pixels and outer pixels */
2206	int j, nrow;
2207
2208	dpsize = pcam->dpixdiameter_cm;
2209	dtsize = dpsize * sqrt(3.) / 2.;
2210	dhsize = dpsize / 2.;
2211
2212	/* Loop over central pixels to generate co-ordinates */
2213
2214	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2215
2216	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2217
2218	pcam->dpixsizefactor[ipixno-1] = 1.;
2219
2220	in = 0;
2221
2222	i = 0;
2223	itot_inside_ring = 0;
2224	iring_no = 0;
2225
2226	/* Calculate the number of pixels out to and including the ring containing pixel number */
2227	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2228
2229	while (itot_inside_ring == 0){
2230
2231	iN = 3(i(i+1)) + 1;
2232
2233	if (ipixno <= iN){
2234	iring_no = i;
2235	itot_inside_ring = iN;
2236	}
2237
2238	i++;
2239	}
2240
2241
2242	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2243
2244	ipix_in_ring = 0;
2245	for (i = 0; i < iring_no; ++i){
2246
2247	ipix_in_ring = ipix_in_ring + 6;
2248	}
2249
2250	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2251	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2252	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2253
2254	isegment = 0;
2255	fsegment_fract = 0.;
2256	if (iring_no > 0) {
2257
2258	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2259
2260	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2261
2262	}
2263
2264	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2265	/* pixel width (flat to flat) dpsize. */
2266
2267	dxfirst_pix = dpsize*iring_no;
2268	dyfirst_pix = 0.;
2269
2270	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2271
2272	ddxseg1 = - dhsize*iring_no;
2273	ddyseg1 = dtsize*iring_no;
2274	ddxseg2 = -dpsize*iring_no;
2275	ddyseg2 = 0.;
2276	ddxseg3 = ddxseg1;
2277	ddyseg3 = -ddyseg1;
2278	ddxseg4 = -ddxseg1;
2279	ddyseg4 = -ddyseg1;
2280	ddxseg5 = -ddxseg2;
2281	ddyseg5 = 0.;
2282	ddxseg6 = -ddxseg1;
2283	ddyseg6 = ddyseg1;
2284
2285	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2286	/* anti-clockwise around the ring by adding the segment to segment vectors. */
2287
2288	switch (isegment) {
2289
2290	case 0:
2291
2292	pcam->dxc[ipixno-1] = 0.;
2293	pcam->dyc[ipixno-1] = 0.;
2294
2295	case 1:
2296	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2297	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2298
2299	break;
2300
2301	case 2:
2302
2303	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2304	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2305
2306	break;
2307
2308	case 3:
2309
2310	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
2311	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
2312
2313	break;
2314
2315	case 4:
2316
2317	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
2318	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
2319
2320	break;
2321
2322	case 5:
2323
2324	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
2325	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2326
2327	break;
2328
2329	case 6:
2330
2331	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2332	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2333
2334	break;
2335
2336	default:
2337
2338	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2339	return(0);
2340
2341	} /* end switch */
2342
2343	} /* end for */
2344
2345	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2346	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2347
2348	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2349
2350	j = pcam->inumcentralpixels;
2351
2352	for(i=0; i<pcam->inumgappixels; i=i+6){
2353	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2354	pcam->dyc[j + i ] = dstarty;
2355	pcam->dpixsizefactor[j + i] = 1.;
2356	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2357	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2358	pcam->dpixsizefactor[j + i + 1] = 1.;
2359	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2360	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2361	pcam->dpixsizefactor[j + i+ 2] = 1.;
2362	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2363	pcam->dyc[j + i + 3] = dstarty;
2364	pcam->dpixsizefactor[j + i+ 3] = 1.;
2365	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2366	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2367	pcam->dpixsizefactor[j + i+ 4] = 1.;
2368	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2369	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2370	pcam->dpixsizefactor[j + i + 5] = 1.;
2371	} /* end for */
2372	} /* end if */
2373
2374	/* generate positions of the outer pixels */
2375
2376	if( pcam->inumbigpixels > 0 ){
2377
2378	j = pcam->inumcentralpixels + pcam->inumgappixels;
2379
2380	for(i=0; i<pcam->inumbigpixels; i++){
2381	pcam->dpixsizefactor[j + i] = 2.;
2382	}
2383
2384	in = 0;
2385
2386	nrow = (int) ceil(dstartx / 2. / dpsize);
2387
2388	while(in < pcam->inumbigpixels){
2389
2390	pcam->dxc[j + in] = dstartx + dpsize;
2391	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2392	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2393	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2394	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2395	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2396	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2397	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2398	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2399	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2400	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2401	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2402	for(i=1; i<nrow; i++){
2403	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2404	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2405	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2406	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2407	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2408	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2409	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2410	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2411	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
2412	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
2413	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
2414	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
2415	}
2416	in = in + 6 * nrow;
2417	dstartx = dstartx + 2. * dpsize;
2418	nrow = nrow + 1;
2419	} /* end while */
2420
2421	} /* end if */
2422
2423	return(pcam->inumpixels);
2424
2425	}
2426	//!@}
2427
2428	//!-----------------------------------------------------------
2429	// @name bpoint_is_in_pix
2430	//
2431	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
2432	//
2433	// @var *pcam structure camera containing all the
2434	// camera information
2435	// @var dx, dy point coordinates in centimeters
2436	// @var ipixnum pixel number (starting at 0)
2437	// @return TRUE if the point is inside the pixel, FALSE otherwise
2438	//
2439	// DP
2440	//
2441	// @date Thu Oct 14 16:59:04 CEST 1999
2442	//------------------------------------------------------------
2443	// @function
2444
2445	//!@{
2446
2447	/****** bpoint_is_in_pix() *************************************/
2448
2449	int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
2450	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
2451	/* the pixel is assumed to be a "closed set" */
2452
2453	double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
2454	double c, xx, yy; /* auxiliary variable */
2455
2456	b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
2457	a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
2458	c = 1./sqrt(3.);
2459
2460	if((ipixnum < 0)\|\|(ipixnum >= pcam->inumpixels)){
2461	fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
2462	fprintf(stderr, "Exiting.\n");
2463	exit(203);
2464	}
2465	xx = dx - pcam->dxc[ipixnum];
2466	yy = dy - pcam->dyc[ipixnum];
2467
2468	if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) \|\|
2469	((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
2470	return(TRUE); /* point is inside */
2471	}
2472	else{
2473	return(FALSE); /* point is outside */
2474	}
2475	}
2476
2477	//!@}
2478
2479	//------------------------------------------------------------
2480	// @name dist_r_P
2481	//
2482	// @desc distance straight line r - point P
2483	//
2484	// @date Sat Jun 27 05:58:56 MET DST 1998
2485	// @function @code
2486	//------------------------------------------------------------
2487	// dist_r_P
2488	//
2489	// distance straight line r - point P
2490	//------------------------------------------------------------
2491
2492	float
2493	dist_r_P(float a, float b, float c,
2494	float l, float m, float n,
2495	float x, float y, float z)
2496	{
2497	return (
2498	sqrt((SQR((a-x)m-(b-y)l) +
2499	SQR((b-y)n-(c-z)m) +
2500	SQR((c-z)l-(a-x)n))/
2501	(SQR(l)+SQR(m)+SQR(n))
2502	)
2503	);
2504	}
2505
2506	//------------------------------------------------------------
2507	// @name check_reflector_file
2508	//
2509	// @desc check if a given reflector file has the right signature
2510	// @desc return TRUE or FALSE
2511	//
2512	// @date Mon Feb 14 16:44:21 CET 2000
2513	// @function @code
2514	//------------------------------------------------------------
2515
2516	int check_reflector_file(FILE *infile){
2517
2518	char Signature[20]; // auxiliary variable
2519	char sign[20]; // auxiliary variable
2520
2521	strcpy(Signature, REFL_SIGNATURE);
2522
2523	strcpy(sign, Signature);
2524
2525	fread( (char *)sign, strlen(Signature), 1, infile);
2526
2527	if (strcmp(sign, Signature) != 0) {
2528	cout << "ERROR: Signature of .rfl file is not correct\n";
2529	cout << '"' << sign << '"' << '\n';
2530	cout << "should be: " << Signature << '\n';
2531	return(FALSE);
2532	}
2533
2534	fread( (char *)sign, 1, 1, infile);
2535
2536	return(TRUE);
2537
2538	}
2539
2540	//------------------------------------------------------------
2541	// @name lin_interpol
2542	//
2543	// @desc interpolate linearly between two points returning the
2544	// @desc y-value of the result
2545	//
2546	// @date Thu Feb 17 11:31:32 CET 2000
2547	// @function @code
2548	//------------------------------------------------------------
2549
2550	float lin_interpol( float x1, float y1, float x2, float y2, float x){
2551
2552	if( (x2 - x1)<=0. ){ // avoid division by zero, return average
2553	cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
2554	return((y1+y2)/2.);
2555	}
2556	else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
2557	return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
2558	}
2559	}
2560
2561
2562	//------------------------------------------------------------
2563	// @name produce_phes
2564	//
2565	// @desc read the photons of an event, pixelize them and simulate QE
2566	// @desc return various statistics and the array of Photoelectrons
2567	//
2568	// @date Mon Feb 14 16:44:21 CET 2000
2569	// @function @code
2570	//------------------------------------------------------------
2571
2572	int produce_phes( FILE *sp, // the input file
2573	struct camera *cam, // the camera layout
2574	float minwl_nm, // the minimum accepted wavelength
2575	float maxwl_nm, // the maximum accepted wavelength
2576	class MTrigger *trigger, // the generated phes
2577	class MFadc *fadc,
2578	int *itotnphe, // total number of produced photoelectrons
2579	float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
2580	int *incph, // total number of cph read
2581	float *tmin_ns, // minimum arrival time of all phes
2582	float *tmax_ns // maximum arrival time of all phes
2583	){
2584
2585	static int i, k, ipixnum;
2586	static float cx, cy, wl, qe, t;
2587	static MCCphoton photon;
2588	static float **qept;
2589	static char flag[SIZE_OF_FLAGS + 1];
2590	static float radius;
2591
2592
2593	// reset variables
2594
2595	for ( i=0; i<cam->inumpixels; ++i ){
2596
2597	nphe[i] = 0.0;
2598
2599	}
2600
2601	*incph = 0;
2602
2603	radius = cam->dxc[cam->inumpixels-1]
2604	+ 1.5cam->dpixdiameter_cmcam->dpixsizefactor[cam->inumpixels-1];
2605
2606	//- - - - - - - - - - - - - - - - - - - - - - - - -
2607	// read photons and "map" them into the pixels
2608	//--------------------------------------------------
2609
2610	// initialize CPhoton
2611
2612	photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
2613
2614	// read the photons data
2615
2616	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2617
2618	// loop over the photons
2619
2620	while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
2621
2622	memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
2623
2624	fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
2625
2626	// increase number of photons
2627
2628	(*incph)++;
2629
2630	// Chceck if photon is inside trigger time range
2631
2632	t = photon.get_t() ;
2633
2634	if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
2635
2636	// read next Photon
2637
2638	fread( flag, SIZE_OF_FLAGS, 1, sp );
2639
2640	// go to beginning of loop, the photon is lost
2641	continue;
2642	}
2643
2644	//
2645	// Pixelization
2646	//
2647
2648	cx = photon.get_x();
2649	cy = photon.get_y();
2650
2651	// get wavelength
2652
2653	wl = photon.get_wl();
2654
2655	// cout << "wl " << wl << " radius " << sqrt(cxcx + cycy) << "\n";
2656
2657	// check if photon has valid wavelength and is inside outermost camera radius
2658
2659	if( (wl > maxwl_nm) \|\| (wl < minwl_nm) \|\| (sqrt(cxcx + cycy) > radius ) ){
2660
2661	// cout << " lost first check\n";
2662
2663	// read next CPhoton
2664
2665	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2666
2667	// go to beginning of loop, the photon is lost
2668	continue;
2669
2670	}
2671
2672	ipixnum = -1;
2673
2674	for(i=0; i<cam->inumpixels; i++){
2675	if( bpoint_is_in_pix( cx, cy, i, cam) ){
2676	ipixnum = i;
2677	break;
2678	}
2679	}
2680
2681	if(ipixnum==-1){// the photon is in none of the pixels
2682
2683	// cout << " lost pixlization\n";
2684
2685	// read next CPhoton
2686
2687	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2688
2689	// go to beginning of loop, the photon is lost
2690	continue;
2691	}
2692
2693	if(ipixnum==0) {// the phton is in the central pixel, which is not used for trigger
2694	// read next CPhoton
2695
2696	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2697
2698	// go to beginning of loop, the photon is lost
2699	continue;
2700	}
2701	//+++
2702	// QE simulation
2703	//---
2704
2705	// set pointer to the QE table of the relevant pixel
2706
2707	qept = (float **)QE[ipixnum];
2708
2709	// check if wl is inside table; outside the table, QE is assumed to be zero
2710
2711
2712	if((wl < qept[0][0]) \|\| (wl > qept[0][pointsQE-1])){
2713
2714	// cout << " lost\n";
2715
2716	// read next Photon
2717
2718	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2719
2720	// go to beginning of loop
2721	continue;
2722
2723	}
2724
2725	// find data point in the QE table (-> k)
2726
2727	k = 1; // start at 1 because the condition was already tested for 0
2728	while (k < pointsQE-1 && qept[0][k] < wl){
2729	k++;
2730	}
2731
2732	// calculate the qe value between 0. and 1.
2733
2734	qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
2735
2736	// if random > quantum efficiency, reject it
2737
2738	if ( (RandomNumber) > qe ) {
2739
2740	// cout << " lost\n";
2741
2742	// read next Photon
2743
2744	fread ( flag, SIZE_OF_FLAGS, 1, sp );
2745
2746	// go to beginning of loop
2747	continue;
2748
2749	}
2750
2751	//+++
2752	// The photon has produced a photo electron
2753	//---
2754
2755	// cout << " accepted\n";
2756
2757	// increment the number of photoelectrons in the relevant pixel
2758
2759	nphe[ipixnum] += 1.0;
2760
2761	// store the new photoelectron
2762
2763	fadc->Fill(ipixnum,(t-tmin_ns) , trigger->FillShow(ipixnum,t-tmin_ns));
2764
2765	*itotnphe += 1;
2766
2767	// read next Photon
2768
2769	fread( flag, SIZE_OF_FLAGS, 1, sp );
2770
2771	} // end while, i.e. found end of event
2772
2773	return(0);
2774
2775	}
2776
2777
2778	//------------------------------------------------------------
2779	// @name produce_nsbrates
2780	//
2781	// @desc read the starfield file, call produce_phes on it in,
2782	// @desc different wavebands, calculate the nsbrates
2783	//
2784	// @date Mon Feb 14 16:44:21 CET 2000
2785	// @function @code
2786	//------------------------------------------------------------
2787
2788	int produce_nsbrates( char *iname, // the starfield input file name
2789	struct camera *cam, // camera layout
2790	float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
2791	// the NSB rates in phe/ns for each pixel
2792	){
2793
2794	int i, j, k, ii; // counters
2795
2796	MTrigger trigger(Trigger_gate_length, Trigger_response_ampl, Trigger_response_fwhm);
2797	MFadc flashadc;
2798
2799	static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
2800	WAVEBANDBOUND2,
2801	WAVEBANDBOUND3,
2802	WAVEBANDBOUND4,
2803	WAVEBANDBOUND5,
2804	WAVEBANDBOUND6 };
2805
2806	FILE *infile; // the input file
2807	fpos_t fileposition; // marker on the input file
2808	static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
2809	static MCEventHeader evth; // the event header
2810	static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
2811	int itnphe; // total number of produced photoelectrons
2812	int itotnphe; // total number of produced photoelectrons after averaging
2813	int incph; // total number of cph read
2814	float tmin_ns; // minimum arrival time of all phes
2815	float tmax_ns; // maximum arrival time of all phes
2816	float integtime_ns; // integration time from the starfield generator
2817
2818	// open input file
2819
2820	log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
2821
2822	infile = fopen( iname, "r" );
2823
2824	// check if the satrfield input file exists
2825
2826	if ( infile == NULL ) {
2827
2828	log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
2829
2830	log( SIGNATURE, "There is not NSB from the Stars\n");
2831
2832	return (0);
2833	}
2834
2835	// get signature, and check it
2836
2837	if(check_reflector_file(infile)==FALSE){
2838	exit(1);
2839	}
2840
2841	// initialize flag
2842
2843	strcpy( cflag, " \0" );
2844
2845	// get flag
2846
2847	fread( cflag, SIZE_OF_FLAGS, 1, infile );
2848
2849	if ( ! feof(infile)){
2850
2851	// reading .rfl file
2852
2853	if(!isA( cflag, FLAG_START_OF_RUN )){
2854	error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
2855	}
2856	else { // found start of run
2857
2858	fread( cflag, SIZE_OF_FLAGS, 1, infile );
2859
2860	if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
2861
2862	// get MCEventHeader
2863
2864	fread( (char*)&evth, evth.mysize(), 1, infile );
2865
2866	integtime_ns = evth.get_energy();
2867
2868	// memorize where we are in the file
2869
2870	if (fgetpos( infile, &fileposition ) != 0){
2871	error( SIGNATURE, "Cannot position in file ...\n");
2872	}
2873
2874	// loop over the wavebands
2875
2876	for(i=0; i<iNUMWAVEBANDS; i++){
2877
2878	// initialize the rate array
2879
2880	for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
2881	rate_phepns[j][i] = 0.;
2882	}
2883
2884	itotnphe = 0;
2885
2886	// read the photons and produce the photoelectrons
2887	// - in order to average over the QE simulation, call the
2888	// production function iNUMNSBPRODCALLS times
2889
2890	for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
2891
2892	// position the file pointer to the beginning of the photons
2893
2894	fsetpos( infile, &fileposition);
2895
2896	// produce photoelectrons
2897
2898	k = produce_phes( infile,
2899	cam,
2900	wl_nm[i],
2901	wl_nm[i+1],
2902	&trigger, // this is a dummy here
2903	&flashadc, // this is a dummy here
2904	&itnphe,
2905	nphe, // we want this!
2906	&incph,
2907	&tmin_ns,
2908	&tmax_ns );
2909
2910	if( k != 0 ){ // non-zero returnvalue means error
2911	cout << "Exiting.\n";
2912	exit(1);
2913	}
2914
2915	for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
2916	rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
2917	}
2918
2919	itotnphe += itnphe;
2920
2921	} // end for(ii=0 ...
2922
2923	fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
2924	wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
2925
2926	} // end for(i=0 ...
2927
2928	}
2929	else{
2930	cout << "Starfield file contains no event.\nExiting.\n";
2931	exit(1);
2932	} // end if( isA ... event
2933	} // end if ( isA ... run
2934	}
2935	else{
2936	cout << "Starfield file contains no run.\nExiting.\n";
2937	exit(1);
2938	}
2939
2940	fclose( infile );
2941
2942	return(0);
2943
2944	}
2945
2946
2947	//------------------------------------------------------------
2948	// @name produce_nsb_phes
2949	//
2950	// @desc produce the photoelectrons from the nsbrates
2951	//
2952	// @date Thu Feb 17 17:10:40 CET 2000
2953	// @function @code
2954	//------------------------------------------------------------
2955
2956	int produce_nsb_phes( float *atmin_ns,
2957	float *atmax_ns,
2958	float theta_rad,
2959	struct camera *cam,
2960	float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
2961	float difnsbr_phepns[iMAXNUMPIX],
2962	float extinction[iNUMWAVEBANDS],
2963	float fnpx[iMAXNUMPIX],
2964	MTrigger *trigger,
2965	MFadc *fadc,
2966	int *inphe,
2967	float base_mv[iMAXNUMPIX]){
2968
2969	float simtime_ns; // NSB simulation time
2970	int i, j, k, ii;
2971	float zenfactor; // correction factor calculated from the extinction
2972	int inumnsbphe; // number of photoelectrons caused by NSB
2973
2974	float t;
2975
2976	ii = *inphe; // avoid dereferencing
2977
2978	// check if the arrival times are set; if not generate them
2979
2980	if(atmin_ns <SIMTIMEOFFSET_NS \|\| atmin_ns > *atmax_ns){
2981	*atmin_ns = 0.;
2982	*atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
2983
2984	}
2985	else{ // extend the event time window by the given offsets
2986
2987	atmin_ns = atmin_ns - SIMTIMEOFFSET_NS;
2988	atmax_ns = atmax_ns + SIMTIMEOFFSET_NS;
2989
2990	simtime_ns = atmax_ns - atmin_ns;
2991
2992	// make sure the simulated time is long enough for the FADC
2993	// simulation and not too long
2994
2995	if(simtime_ns< TOTAL_TRIGGER_TIME){
2996	atmin_ns = atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
2997	atmax_ns = atmin_ns + TOTAL_TRIGGER_TIME;
2998	simtime_ns = TOTAL_TRIGGER_TIME;
2999	}
3000
3001	if(simtime_ns> TOTAL_TRIGGER_TIME){
3002	atmax_ns =atmin_ns + TOTAL_TRIGGER_TIME;
3003	simtime_ns = TOTAL_TRIGGER_TIME;
3004	}
3005
3006	}
3007
3008	// initialize baselines
3009
3010	for(i=0; i<cam->inumpixels; i++){
3011	base_mv[i] = 0.;
3012	}
3013
3014	// calculate baselines and generate phes
3015
3016	for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
3017
3018	// calculate the effect of the atmospheric extinction
3019
3020	zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
3021
3022	for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
3023
3024	inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
3025	* simtime_ns );
3026
3027	base_mv[j] += inumnsbphe;
3028
3029	// randomize
3030
3031	inumnsbphe = ignpoi( inumnsbphe );
3032
3033	// create the photoelectrons
3034
3035	for(k=0; k < inumnsbphe; k++){
3036
3037	t=(RandomNumber * simtime_ns);
3038
3039	(fadc).Fill(j,t ,(trigger).FillNSB(j,t));
3040
3041	ii++; // increment total number of photoelectons
3042
3043	fnpx[j] += 1.; // increment number of photoelectrons in this pixel
3044
3045	}
3046
3047	} // end for(j=0 ...
3048	} // end for(i=0 ...
3049
3050	// finish baseline calculation
3051
3052	for(i=0; i<cam->inumpixels; i++){
3053	base_mv[i] = RESPONSE_FWHM RESPONSE_AMPLITUDE / simtime_ns;
3054	}
3055
3056	*inphe = ii; // update the pointer
3057
3058	return(0);
3059	}
3060
3061
3062	// @endcode
3063
3064
3065	//=------------------------------------------------------------
3066	//!@subsection Log of this file.
3067
3068	//!@{
3069	//
3070	// $Log: not supported by cvs2svn $
3071	// Revision 1.11 2000/09/21 11:47:33 harald
3072	// Oscar found some smaller errors in the calculation of the pixel shape and
3073	// corrected it.
3074	//
3075	// Revision 1.10 2000/07/04 14:10:20 MagicSol
3076	// Some changes have been done in the root output file. The RawEvt tree is only
3077	// stored in the single trigger mode.
3078	// The trigger input parameters are also given by the general input card.
3079	// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
3080	// each one can be studied seperately.
3081	//
3082	// Revision 1.9 2000/06/13 13:25:24 blanch
3083	// The multiple arrays have been replaced, since they do not work
3084	// in alpha machines. Now we are using pointers and the command new
3085	// to allocate memory.
3086	//
3087	// Revision 1.8 2000/05/11 13:57:27 blanch
3088	// The option to loop over several trigger configurations has been included.
3089	// Some non-sense with triggertime range has been solved.
3090	// Montecarlo information and ADC counts are saved in a root file.
3091	// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
3092	//
3093	// Revision 1.7 2000/03/24 18:10:46 blanch
3094	// A first FADC simulation and a trigger simulation are already implemented.
3095	// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
3096	// 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.
3097	//
3098	// Revision 1.6 2000/03/20 18:35:11 blanch
3099	// 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.
3100	//
3101	// Revision 1.5 2000/02/18 17:40:35 petry
3102	// This version includes drastic changes compared to camera.cxx 1.4.
3103	// It is not yet finished and not immediately useful because the
3104	// trigger simulation is not yet re-implemented. I had to take it
3105	// out together with some other stuff in order to tidy the whole
3106	// program up. This is not meant as an insult to anyone. I needed
3107	// to do this in order to be able to work on it.
3108	//
3109	// This version has been put in the repository in order to be
3110	// able to share the further development with others.
3111	//
3112	// If you need something working, wait or take an earlier one.
3113	// See file README.
3114	//
3115	// Revision 1.4 2000/01/25 08:36:23 petry
3116	// The pixelization in previous versions was buggy.
3117	// This is the first version with a correct pixelization.
3118	//
3119	// Revision 1.3 2000/01/20 18:22:17 petry
3120	// Found little bug which makes camera crash if it finds a photon
3121	// of invalid wavelength. This bug is now fixed and the range
3122	// of valid wavelengths extended to 290 - 800 nm.
3123	// This is in preparation for the NSB simulation to come.
3124	// Dirk
3125	//
3126	// Revision 1.2 1999/11/19 08:40:42 harald
3127	// Now it is possible to compile the camera programm under osf1.
3128	//
3129	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
3130	// This the starting point for CVS controlled further developments of the
3131	// camera program. The program was originally written by Jose Carlos.
3132	// But here you can find a "rootified" version to the program. This means
3133	// that there is no hbook stuff in it now. Also the output of the
3134	// program changed to the MagicRawDataFormat.
3135	//
3136	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
3137	//
3138	// Revision 1.3 1999/10/22 15:01:28 petry
3139	// version sent to H.K. and N.M. on Fri Oct 22 1999
3140	//
3141	// Revision 1.2 1999/10/22 09:44:23 petry
3142	// first synthesized version which compiles and runs without crashing;
3143	//
3144	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
3145	// first synthesised version
3146	//
3147	// Revision 1.13 1999/03/15 14:59:05 gonzalez
3148	// camera-1_1
3149	//
3150	// Revision 1.12 1999/03/02 09:56:10 gonzalez
3151	// * empty log message *
3152	//
3153	//
3154	//!@}
3155
3156	//=EOF

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

Download in other formats: