Context Navigation

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

Visit:

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

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

Download in other formats: