Context Navigation

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

Visit:

Last change on this file since 1198 was 1195, checked in by blanch, 23 years ago
The option of adding noise to all pixels or to not adding the noise has been added. We removed the pixels larger than 577. When there were more than one trigger in one shower, the pixel number was increasing. Now it is flagged by the variable MMcTrig::fFirstLvlTrig.
File size: 98.5 KB

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

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

Download in other formats: