Context Navigation

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

Visit:

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

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

Download in other formats: