Context Navigation

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

Visit:

Last change on this file since 2427 was 2427, checked in by blanch, 21 years ago
- The screen output information has been improved to prevent some non-desired running conditions just looking at the output screen. - One MMcEvt for each Telscope is stored in the output file. - 500 empty events are simualted to get a more precise estimation of the pedestal Sigma for each pixel.
File size: 127.9 KB

Line
1	//!/////////////////////////////////////////////////////////////////////
2	//
3	// camera
4	//
5	// @file camera.cxx
6	// @title Camera simulation
7	// @subtitle Code for the simulation of the camera phase
8	// @desc Code for the simulation of the camera of CT1 and MAGIC
9	// @author J C Gonzalez
10	// @email gonzalez@mppmu.mpg.de
11	// @date Thu May 7 16:24:22 1998
12	//
13	//----------------------------------------------------------------------
14	//
15	// Created: Thu May 7 16:24:22 1998
16	// Author: Jose Carlos Gonzalez
17	// Purpose: Program for reflector simulation
18	// Notes: See files README for details
19	//
20	//----------------------------------------------------------------------
21	//
22	// $RCSfile: camera.cxx,v $
23	// $Revision: 1.65 $
24	// $Author: blanch $
25	// $Date: 2003-10-26 19:43:00 $
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	#include "TObjArray.h"
61
62	#include "MTrigger.hxx"
63	#include "MFadc.hxx"
64	#include "MLons.hxx"
65
66	#include "MRawRunHeader.h"
67	#include "MRawEvtData.h"
68	#include "MRawEvtHeader.h"
69	#include "MRawCrateData.h"
70	#include "MMcEvt.hxx"
71	#include "MMcTrig.hxx"
72	#include "MMcRunHeader.hxx"
73	#include "MMcConfigRunHeader.h"
74	#include "MMcCorsikaRunHeader.h"
75	#include "MMcTrigHeader.hxx"
76	#include "MMcFadcHeader.hxx"
77	#include "MGeomCamMagic.h"
78	#include "MGeomCamMagic919.h"
79	#include "MGeomCamMagicHG.h"
80	#include "MGeomCamECO1000.h"
81	#include "MGeomCamECO1000HG.h"
82	#include "MGeomCamCT1.h"
83	#include "MGeomCamCT1Daniel.h"
84	#include "MGeomPix.h"
85
86	/*!@"
87
88	All the defines are located in the file \|camera.h\|.
89
90	@"*/
91
92	#include "camera.h"
93
94	//!@}
95
96	/*!@"
97
98	The following set of flags are used in time of compilation. They do
99	not affect directly the behaviour of the program at run-time
100	(though, of course, if you disconnected the option for
101	implementation of the Trigger logic, you will not be able to use any
102	trigger at all. The 'default' values mean default in the sense of
103	what you got from the server when you obtained this program.
104
105	@"*/
106
107	//!@{
108
109	// flag for debugging (default: OFF )
110	#define __DEBUG__
111	#undef __DEBUG__
112
113
114	//!@}
115
116	//=-----------------------------------------------------------
117	//!@subsection Definition of global variables.
118
119	/*!@"
120
121	Now we define some global variables with data about the telescope,
122	such as "focal distance", number of pixels/mirrors,
123	"size of the camera", and so on.
124
125	@"*/
126
127	/*!@"
128
129	Depending on the telescope we are using (CT1 or MAGIC), the
130	information stored in the definition file is different.
131	The variable \|ct_Type\| has the value 0 when we use
132	CT1, and 1 when we use MAGIC.
133
134	@"*/
135
136	/*!@"
137
138	And this is the information about the whole telescope.
139
140	@"*/
141
142	//!@{
143
144	// parameters of the CT (from the CT definition file)
145
146	//@: Number of pixels
147	static int ct_NPixels;
148
149	//@: Number of CT
150	static int ct_Number;
151
152	//@: Camera geometries
153	static int ct_Geometry;
154
155	//@: list of showers to be skipped
156	static int *Skip;
157
158	//@: number of showers to be skipped
159	static int nSkip=0;
160
161	//@: flag: TRUE: data come from STDIN; FALSE: from file
162	static int Data_From_STDIN = FALSE;
163
164	//@: flag: TRUE: write all images to output; FALSE: only triggered showers
165	static int Write_All_Images = FALSE;
166
167	static int Write_McEvt = TRUE;
168	static int Write_McTrig = TRUE;
169	static int Write_McFADC = TRUE;
170	static int Write_RawEvt = FALSE;
171
172	//@: flag: TRUE: selection on the energy
173	static int Select_Energy = TRUE;
174
175	//@: Lower edge of the selected energy range (in GeV)
176	static float Select_Energy_le = 0.0;
177
178	//@: Upper edge of the selected energy range (in GeV)
179	static float Select_Energy_ue = 100000.0;
180
181	//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
182	static int FADC_Scan = FALSE;
183
184	//@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't
185	static int Trigger_Scan = FALSE;
186
187	//@: flag: TRUE: loop trigger analysis over several thresholds, multiplicities and topologies; FALSE: a single trigger configuration
188	static int Trigger_Loop = FALSE;
189
190	//@: flag: TRUE: Different threshold for each pixel ; FALSE: same threshold for all pixels
191	static int Individual_Thres_Pixel = FALSE;
192
193	//@: Properties of the trigger
194	static float Trigger_gate_length = 6.0;
195	static float Trigger_response_ampl = 1.0;
196	static float Trigger_response_fwhm = 2.0;
197	static float Trigger_overlaping_time= 0.25;
198	static float Trigger_noise= 0.3;
199
200	//@: Properties of the FADC
201	static float FADC_response_ampl = MFADC_RESPONSE_AMPLITUDE;
202	static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
203	static float FADC_resp_ampl_out = MFADC_RESPONSE_AMPLITUDE;
204	static float FADC_resp_fwhm_out = MFADC_RESPONSE_FWHM;
205	static float FADC_noise = 2.0;
206	static float DIGITAL_noise = 2.0;
207	static float FADC_high2low = HIGH2LOWGAIN;
208
209	//@: Trigger conditions for a single trigger mode
210	static float **qThreshold;
211	static int Trigger_multiplicity[100];
212	static int Trigger_topology[100];
213
214	//@: Upper and lower edges of the trigger loop
215	static float Trigger_loop_lthres = 2.0;
216	static float Trigger_loop_uthres = 10.0;
217	static float Trigger_loop_sthres = 1.0;
218	static int Trigger_loop_lmult = 2;
219	static int Trigger_loop_umult = 10;
220	static int Trigger_loop_ltop = 0;
221	static int Trigger_loop_utop = 2;
222
223	//!@}
224
225	/*!@"
226
227	The following double-pointer is a 2-dimensional table with information
228	about each pixel. The routine read_pixels will generate
229	the information for filling it using igen_pixel_coordinates().
230
231	@"*/
232
233	//!@{
234	// Pointer to a tables/Arrays with information about the pixels
235	// and data stored on them with information about the pixels
236
237
238	//@: contents of the pixels (ph.e.)
239	static float *fnpix;
240
241
242	//!@}
243
244	/*!@"
245
246	The following double-pointer is a 2-dimensional table with the
247	Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
248	of the wavelength @$\lambda@$. The routine \|read_pixels()\| will read
249	also this information from the file \|qe.dat\|.
250
251	@"*/
252
253	//!@{
254	// Pointer to a table with QE, number of datapoints, and wavelengths
255
256	//@: table of QE
257	static float ****QE;
258
259	//@: number of datapoints for the QE curve
260	static int pointsQE[100];
261
262	//@: table of QE
263	static float *QElambda;
264
265	//@: table of lightguide = WC;
266	static float **WC;
267
268	//@: number of datapoints for the WC curve
269	static int pointsWC;
270
271	//!@}
272
273	/*!@"
274
275	The following double-pointer is a 2-dimensional table with information
276	about each mirror in the dish. The routine \|read_ct_file()\| will read
277	this information from the CT definition file.
278
279	@"*/
280
281	/*!@"
282
283	We define a table into where random numbers will be stored.
284	The routines used for random number generation are provided by
285	\|RANLIB\| (taken from NETLIB, \|www.netlib.org\|), and by
286	the routine \|double drand48(void)\| (prototype defined in
287	\|stdlib.h\|) through the macro \|RandomNumber\| defined in
288	\|camera.h\|.
289
290	@"*/
291
292
293	//!@}
294
295	extern char FileName[];
296	// switch on starfield rotation
297	static int Starfield_rotate = TRUE;
298
299	//=-----------------------------------------------------------
300	// @subsection Main program.
301
302	//!@{
303
304	//++++++++++++++++++++++++++++++++++++++++
305	// MAIN PROGRAM
306	//----------------------------------------
307
308	int main(int argc, char **argv)
309	{
310
311	//!@' @#### Definition of variables.
312	//@'
313
314	char inname_CT[100][256]; //@< array of names for each CT input file
315	char starfieldname[256]; //@< starfield input file name
316	char qe_filename[256]; //@< qe input file name
317
318	char datname[256]; //@< data (ASCII) output file name
319
320	char rootname[256] ; //@< ROOT file name
321	char rootname_loop[256] ; //@< ROOT file name
322
323	char parname[256]; //@< parameters file name
324
325	char nsbpathname[256]; //@< directory with the dataabse for th NSB
326	char nsbpath_outer[256]; //@< directory with the dataabse for outer pixels
327
328	char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
329	char flag_new[SIZE_OF_FLAGS + 1]; //@< New flag for the run header in the .rfl file
330	char GeometryName[100][50]; //@< Name of MGeomCam for each CT
331	int GeometryCamera[100]; //@< Identification of MGeomCam for each CT
332	int TriggerPixels[100]; //@< Numeber of pixels in the trigger region for each CT
333
334	int reflector_file_version=0; //@< vrsion of he reflector file
335
336	FILE *inputfile[100]; //@< stream for the input file
337
338	ofstream datafile; //@< stream for the data file
339
340	MCRunHeader mcrunh; //@< Run Header class (MC)
341	MCEventHeader mcevth[100]; //@< Event Header class (MC)
342	MCEventHeader_2 mcevth_2[100]; //@< Event Header class (MC) for reflector > 0.6
343	MCCphoton cphoton; //@< Cherenkov Photon class (MC)
344
345	int inumphe; //@< number of photoelectrons in an event from showers
346	int inumphe_CT[100]; //@< number of photoelectrons in an event from showers
347	float inumphensb[100]; //@< number of photoelectrons in an event from nsb
348
349	float arrtmin_ns; //@ arrival time of the first photoelectron
350	float arrtmax_ns; //@ arrival time of the last photoelectron
351
352	float thetaCT[100], phiCT[100]; //@< parameters of a given shower
353	float thetashw, phishw; //@< parameters of a given shower
354	float coreX, coreY; //@< core position
355	float impactD[100]; //@< impact parameter
356	float l1, m1, n1; //@< auxiliary variables
357	float l2, m2, n2; //@< auxiliary variables
358	float num, den; //@< auxiliary variables
359	float factorqe_NSB[100]; //@< factor on the NSB depending of QE
360
361	int nshow=0; //@< partial number of shower in a given run
362	int ntshow=0; //@< total number of showers
363	int ncph[100]; //@< partial number of photons in a given run
364	int ntcph[100]; //@< total number of photons
365
366	int ibr, j, k; //@< simple counters
367
368	int addElecNoise; //@< Will we add ElecNoise?
369
370	float riseDiskThres;
371	float secureDiskThres;
372
373	int simulateNSB; //@< Will we simulate NSB?
374	int nphe2NSB; //@< From how many phe will we simulate NSB?
375	float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel)
376	float diffnsb_phepns[100][iMAXNUMPIX]; //@< diffuse NSB values for each pixel
377	//@< derived from meanNSB
378	float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS]; //@< non-diffuse nsb
379	//@< photoelectron rates
380	float nsb_phepns[100][iMAXNUMPIX]; //@< NSB values for each pixel
381	float nsb_phepns_rotated[100][iMAXNUMPIX]; //@< NSB values for each pixel after rotation
382
383	Float_t nsb_trigresp[TRIGGER_TIME_SLICES]; //@< array to write the trigger
384	//@< response from the database
385	Float_t nsb_fadcresp[(Int_t) SLICES_MFADC]; //@< array to write the fadc
386	//@< response from the database
387	Byte_t trigger_map[((Int_t)(CAMERA_PIXELS/8))+1]; //@< Pixels on when the
388	//@< camera triggers
389	Float_t fadc_elecnoise[CAMERA_PIXELS]; //@< Electronic noise for each pixel
390	Float_t fadc_diginoise[CAMERA_PIXELS]; //@< Digital noise for each pixel
391	Float_t fadc_pedestals[CAMERA_PIXELS]; //@< array for fadc pedestals values
392	Float_t fadc_sigma[CAMERA_PIXELS]; //@< array for fadc pedestals sigma
393	Float_t fadc_sigma_low[CAMERA_PIXELS]; //@< array for fadc pedestals sigma
394
395	float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
396	EXTWAVEBAND1,
397	EXTWAVEBAND2,
398	EXTWAVEBAND3,
399	EXTWAVEBAND4,
400	EXTWAVEBAND5
401	};
402	float zenfactor=1.0; // correction factor calculated from the extinction
403
404	int nstoredevents = 0;
405	int flagstoring = 0;
406
407	int ntrigger[100]; //@< number of triggers in the whole file
408	int btrigger = 0; //@< trigger flag
409	int ithrescount; //@< counter for loop over threshold trigger
410	float fthrescount; //@< value for loop over threshold trigger
411	int imulticount; //@< counter for loop over multiplicity trigger
412	int itopocount; //@< counter for loop over topology trigger
413	int isorttopo[3]; //@< sorting the topologies
414	int icontrigger; //@< number of trigger conditions to be analised
415
416	float fpixelthres[CAMERA_PIXELS]; //@< Threshold values
417
418	TArrayC *fadcValues; //@< the analog Fadc High gain signal for pixels
419	TArrayC *fadcValuesLow; //@< the analog Fadc Low gain signal for pixels
420
421	int still_in_loop = FALSE;
422
423	//@< Variables to fill the McRunHeader
424	Int_t sfRaH = 5;
425	Int_t sfRaM = 34;
426	Int_t sfRaS = 32;
427	Int_t sfDeD = 22;
428	Int_t sfDeM = 00;
429	Int_t sfDeS = 52;
430	Float_t shthetamax = 0.0;
431	Float_t shthetamin = 0.0;
432	Float_t shphimax = 0.0;
433	Float_t shphimin = 0.0;
434	Float_t telestheta = 0.0;
435	Float_t telesphi = 0.0;
436	Float_t sofftheta = 0.0;
437	Float_t soffphi = 0.0;
438	UInt_t corsika = 5200 ;
439	Float_t maxpimpact = 0.0;
440
441	// Star Field Rotation variables
442	Float_t zenith = 0.0;
443	Float_t azimutal = 90.0;
444	Float_t rho , C3 , C2 , C1;
445
446	//!@' @#### Definition of variables for \|getopt()\|.
447	//@'
448
449	int ch, errflg = 0; //@< used by getopt
450
451	//!@' @#### Initialisation of some variables
452	//@'
453
454	qThreshold = new float *[100];
455	for(int i=0;i<100;i++)
456	qThreshold[i] = new float [CAMERA_PIXELS];
457
458	for(int i=0;i<iMAXNUMPIX;i++){
459	for(int ict=0;ict<100;ict++){
460	nsb_phepns[ict][i]=0;
461	nsb_phepns_rotated[ict][i]=0;
462	}
463	for(j=0;j<iNUMWAVEBANDS;j++)
464	nsbrate_phepns[i][j]=0.0; //@< Starlight rates initialised at 0.0
465	}
466	for(int i=0;i<100;i++)
467	ntcph[i]=0;
468	/*!@'
469
470	@#### Beginning of the program.
471
472	We start with the main program. First we (could) make some
473	presentation, and follows the reading of the parameters file (now
474	from the \|stdin\|), the reading of the CT parameters file, and the
475	creation of the output file, where the processed data will be
476	stored.
477
478	*/
479
480	//++
481	// START
482	//--
483
484	// make unbuffered output
485
486	cout.setf ( ios::stdio );
487
488	// parse command line options (see reflector.h)
489
490	parname[0] = '\0';
491
492	optarg = NULL;
493	while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
494	switch (ch) {
495	case 'f':
496	strcpy(parname, optarg);
497	break;
498	case 'h':
499	usage();
500	break;
501	default :
502	errflg++;
503	}
504
505	// show help if error
506
507	if ( errflg>0 )
508	usage();
509
510	// make some sort of presentation
511
512	present();
513
514	// read parameters file
515
516	if ( strlen(parname) < 1 )
517	readparam(NULL);
518	else
519	readparam(parname);
520
521	// read data from file or from STDIN?
522
523	Data_From_STDIN = get_data_from_stdin();
524
525	// get number of telescopes and camera geometries
526
527	ct_Number=get_ct_number();
528	ct_Geometry=get_ct_geometry();
529
530	for(int ict=0;ict<ct_Number;ict++){
531	ntrigger[ict]=0;
532	if(ct_Geometry>=ict*10){
533	GeometryCamera[ict]=int(ct_Geometry/pow(10.0,ict))%10;
534	}
535	else
536	GeometryCamera[ict]=0;
537	}
538
539	// structure holding the camera definition
540
541	TObjArray camgeom;
542
543	for (int i=0; i<ct_Number;i++){
544	switch(GeometryCamera[i]){
545	case 1: camgeom[i]=new MGeomCamMagic;
546	strcpy(GeometryName[i],"MGeomCamMagic");
547	TriggerPixels[i]=TRIGGER_PIXELS_1;
548	break;
549	case 2: camgeom[i]=new MGeomCamMagic919;
550	strcpy(GeometryName[i],"MGeomCamMagic919");
551	TriggerPixels[i]=TRIGGER_PIXELS_2;
552	break;
553	case 3: camgeom[i]=new MGeomCamMagicHG;
554	strcpy(GeometryName[i],"MGeomCamMagicHG");
555	TriggerPixels[i]=TRIGGER_PIXELS_3;
556	break;
557	case 5: camgeom[i]=new MGeomCamECO1000;
558	strcpy(GeometryName[i],"MGeomCamECO1000");
559	TriggerPixels[i]=TRIGGER_PIXELS_5;
560	break;
561	case 6: camgeom[i]=new MGeomCamECO1000HG;
562	strcpy(GeometryName[i],"MGeomCamECO1000HG");
563	TriggerPixels[i]=TRIGGER_PIXELS_6;
564	break;
565	case 8: camgeom[i]=new MGeomCamCT1;
566	strcpy(GeometryName[i],"MGeomCamCT1");
567	TriggerPixels[i]=TRIGGER_PIXELS_8;
568	break;
569	case 9: camgeom[i]=new MGeomCamCT1Daniel;
570	strcpy(GeometryName[i],"MGeomCamCT1Daniel");
571	TriggerPixels[i]=TRIGGER_PIXELS_9;
572	break;
573	default:
574	error( SIGNATURE, "Camera geometry %i is not defined. Stoping camera program ...", GeometryCamera[i]);
575	exit(1);
576	break;
577	}
578	}
579
580	ct_NPixels=0;
581
582	for(int ict=0;ict<ct_Number;ict++)
583	ct_NPixels=(((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels()>(UInt_t)ct_NPixels)?
584	(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels():
585	ct_NPixels;
586
587	// read WC and QE files
588
589	// FIX ME !
590	// Currently the PMT N of any camera with same average QE has the same QE.
591
592	QE = new float *** [ct_Number];
593	for(int ict=0;ict<ct_Number;ict++){
594	strcpy( qe_filename, get_qe_filename(ict));
595	read_QE(qe_filename,ict);
596	}
597
598	read_WC();
599
600	// write all images, even those without trigger?
601
602	Write_All_Images = get_write_all_events();
603
604	Write_McEvt = get_write_McEvt() ;
605	Write_McTrig = get_write_McTrig() ;
606	Write_McFADC = get_write_McFadc() ;
607	Write_RawEvt = get_write_RawEvt() ;
608
609	FADC_Scan = get_FADC_Scan();
610	Trigger_Scan = get_Trigger_Scan();
611
612	Individual_Thres_Pixel = get_indi_thres_pixel();
613
614	get_secure_threhold(&riseDiskThres,&secureDiskThres);
615
616	get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm,
617	&FADC_resp_ampl_out, &FADC_resp_fwhm_out);
618	FADC_high2low=get_High_to_Low();
619
620	// FIXME --- tirgger properties may depend on the Camrea geometry
621
622	get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);
623
624	Trigger_Loop = get_Trigger_Loop(&Trigger_loop_lthres, &Trigger_loop_uthres, &Trigger_loop_sthres, &Trigger_loop_lmult, &Trigger_loop_umult, &Trigger_loop_ltop, &Trigger_loop_utop);
625
626	icontrigger =((int)((Trigger_loop_uthres-Trigger_loop_lthres)
627	/Trigger_loop_sthres)+1)*
628	(Trigger_loop_umult-Trigger_loop_lmult+1)*
629	(Trigger_loop_utop-Trigger_loop_ltop+1);
630
631	// Trigger loop operation is not implemented for Multi telscopes
632
633	if ( Trigger_Loop && ct_Number > 1 ){
634	cout<<"ERROR : camera::main : Trigger loop option is not";
635	cout<<" implemented for Multi Telescopes option. "<<icontrigger<<
636	" "<<ct_Number<<endl;
637	exit(1);
638	}
639
640	if (!Trigger_Loop){
641	get_Trigger_Single (qThreshold,
642	&Trigger_multiplicity[0],
643	&Trigger_topology[0]);
644	icontrigger=1;
645	}
646	else
647	get_threshold(qThreshold[0]);
648
649	// get filenames
650
651	for(int ict=0;ict<ct_Number;ict++)
652	strcpy( inname_CT[ict], get_input_filename(ict) );
653
654	strcpy( starfieldname, get_starfield_filename() );
655	strcpy( datname, get_data_filename() );
656	strcpy( rootname, get_root_filename() );
657	strcpy( rootname_loop, get_loop_filename() );
658	strcpy( nsbpathname, get_nsb_directory() );
659	strcpy( nsbpath_outer, get_nsb_directory_outer() );
660
661	// get different parameters of the simulation
662
663	addElecNoise = add_elec_noise(&FADC_noise, &DIGITAL_noise, &Trigger_noise);
664	simulateNSB = get_nsb( &meanNSB, &nphe2NSB );
665
666
667	// get selections on the parameters
668
669	Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
670
671	//Parameters for starfield rotation
672
673	Starfield_rotate = get_starfield_rotate();
674
675	// log filenames information
676	for(Int_t ict=0;ict<ct_Number;ict++){
677
678	log(SIGNATURE,"\t%s : %i\n\t%20s:\t%s\n",
679	"------- TELESCOPE ",ict+1,
680	"Geometry ", GeometryName[ict]);
681	strcpy( qe_filename, get_qe_filename(ict));
682
683	// Look to factor for NSB respect to emi_coat PMTs
684	if(!strcmp(qe_filename,"../Data/qe-emi-coat.RFL.dat"))
685	factorqe_NSB[ict]=EMICOAT_NSB;
686	else if(!strcmp(qe_filename,"../Data/qe-emi.RFL.dat"))
687	factorqe_NSB[ict]=EMI_NSB;
688	else if(!strcmp(qe_filename, "../Data/qe-intevac_hpd.RFL.dat"))
689	factorqe_NSB[ict]=HPD_NSB;
690	else{
691	log(SIGNATURE,"WARNING : the factor for the diffuse NSB for this QE file (%s) is not known. The number of phe(ns stated in the inputcard will be used.\n",qe_filename);
692	factorqe_NSB[ict]=1.0;
693	}
694	log(SIGNATURE,
695	"%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s\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",
696	"Filenames",
697	"In", inname_CT[ict],
698	"Stars", starfieldname,
699	"NSB database","(inner pixels)", nsbpathname,
700	"(outer pixels)", nsbpath_outer,
701	"QE", qe_filename,
702	"Data", datname,
703	"ROOT", rootname
704	);
705
706	// log parameters information
707
708	log(SIGNATURE,
709	"%s:\n\t%20s: %f: %s\n\t%20s: %f\n\t%20s: %f\n",
710	"Parameters",
711	"NSB (phes/ ns 0.10.1 deg^2 239 m^2)", meanNSBfactorqe_NSB[ict], ONoff(simulateNSB),
712	"Pedestals = ", get_FADC_pedestal(),
713	"High to Low gain = ", FADC_high2low
714	);
715
716	// log Trigger information
717
718	if (Trigger_Loop) {
719	log(SIGNATURE,
720	"%s:\n\t%20s: from %5.2f to %5.2f by %5.2f step\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s\n",
721	"Trigger Loop mode",
722	"Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres,
723	"Multiplicity",Trigger_loop_lmult,Trigger_loop_umult,
724	"Topology",Trigger_loop_ltop,Trigger_loop_utop,
725	rootname_loop);
726	}
727	else if (Individual_Thres_Pixel == FALSE){
728	log(SIGNATURE,
729	"%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n",
730	"Single Trigger mode",
731	"Threshold",qThreshold[ict][0],
732	"Multiplicity",Trigger_multiplicity[ict],
733	"Topology",Trigger_topology[ict]);
734	}
735	else{
736	log(SIGNATURE,
737	"%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n",
738	"Single Trigger mode",
739	"Threshold","Different for each pixel",
740	"Multiplicity",Trigger_multiplicity[ict],
741	"Topology",Trigger_topology[ict]);
742	}
743	log(SIGNATURE,"\t%s\n",
744	"END TELESCOPE --------");
745	}
746	// log flags information
747
748	log(SIGNATURE,
749	"%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %3.2f(%s) %3.2f+%3.2f(%s) %s\n",
750	"Flags",
751	"Data_From_STDIN", ONoff(Data_From_STDIN),
752	"Write_All_Events", ONoff(Write_All_Images),
753	"Rotate Starfield", ONoff(Starfield_rotate),
754	"Electronic Noise", Trigger_noise,"trigger",
755	FADC_noise,DIGITAL_noise,"fadc",ONoff(addElecNoise)
756	);
757
758	// log flags information
759
760	log(SIGNATURE,
761	"%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
762	"Root output",
763	"Write_McEvt", ONoff(Write_McEvt),
764	"Write_McTrig", ONoff(Write_McTrig),
765	"Write_McFADC", ONoff(Write_McFADC),
766	"Write_RawEvt", ONoff(Write_RawEvt));
767
768	// log selections
769
770	log(SIGNATURE,
771	"%s:\n\t%20s: %s (%f:%f)\n",
772	"Selections:",
773	"Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
774
775	// Definition and initialization of array to save trigger statistics
776
777	int ***ntriggerloop;
778
779	ntriggerloop= new int ** [(int)((Trigger_loop_uthres-Trigger_loop_lthres)
780	/Trigger_loop_sthres)];
781	for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;fthrescount+=Trigger_loop_sthres, ithrescount++){
782	ntriggerloop[ithrescount]= new int * [Trigger_loop_umult-Trigger_loop_lmult+1];
783	for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
784	ntriggerloop[ithrescount][imulticount]= new int [Trigger_loop_utop-Trigger_loop_ltop+1];
785	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
786	ntriggerloop[ithrescount][imulticount][itopocount]=0;
787	}
788	}
789	}
790
791	// We should be careful that topologies are sort from
792	// the less to the more restrictive one.
793
794	if (Trigger_loop_utop==Trigger_loop_ltop)
795	for(int is=0; is<3;is++)
796	isorttopo[is]=is;
797	else {
798	isorttopo[0]=1;
799	isorttopo[1]=0;
800	isorttopo[2]=2;
801	}
802
803	// get list of showers to evt. skip
804
805	nSkip = get_nskip_showers();
806
807	if (nSkip > 0) {
808	Skip = new int[ nSkip ];
809	get_skip_showers( Skip );
810
811	log(SIGNATURE, "There are some showers to skip:\n");
812	for (int i=0; i<nSkip; ++i)
813	log(SIGNATURE, "\tshower # %d\n", Skip[i]);
814	}
815
816
817	// read parameters from the ct.def file
818
819	//read_ct_file();
820
821	Int_t Lev0, Lev1;
822	Int_t Lev0MT[100], Lev1MT[100];
823
824	fadcValues = new TArrayC(FADC_SLICES);
825	fadcValuesLow = new TArrayC(FADC_SLICES);
826
827	// number of pixels for parameters
828
829	// Switched off writing TObject
830
831	MArray::Class()->IgnoreTObjectStreamer();
832	MParContainer::Class()->IgnoreTObjectStreamer();
833
834	// initialise ROOT
835
836	TROOT simple("simple", "MAGIC Telescope Monte Carlo");
837
838	// initialise instance of Trigger and FADC classes
839
840	MTrigger **Trigger_CT;
841	Trigger_CT = new MTrigger *[ct_Number];
842
843	for (int i=0; i<ct_Number;i++){
844	Trigger_CT[i] = new MTrigger(TriggerPixels[i],
845	((MGeomCam*)(camgeom.UncheckedAt(i))),
846	Trigger_gate_length,
847	Trigger_overlaping_time,
848	Trigger_response_ampl,
849	Trigger_response_fwhm); //@< A instance of the Class MTrigger
850	Trigger_CT[i]->SetSeed(UInt_t(i+get_seeds(0)));
851	}
852
853	for (int ict=0; ict<ct_Number;ict++){
854	// Generate databse for trigger electronic noise
855	if (addElecNoise)
856	Trigger_CT[ict]->SetElecNoise(Trigger_noise) ;
857
858	// Set Right Discriminator threshold, taking into account trigger pixels
859	Trigger_CT[ict]->CheckThreshold(&qThreshold[ict][0],GeometryCamera[ict]);
860	}
861	// Set flag in pixel 0 (not used for trigger) that indicates if secure pixel
862	// is active: secureDiskThres*10000+riseDiskThres
863
864	for (int ict=0; ict<ct_Number;ict++){
865	if(riseDiskThres>0.0)
866	qThreshold[ict][0]=((UInt_t)secureDiskThres100)100+riseDiskThres;
867	else
868	qThreshold[ict][0]=0.0;
869	}
870	// Initialise McTrig information class if we want to save trigger informtion
871
872	MMcTrig **McTrig = NULL;
873	MMcTrigHeader **HeaderTrig = NULL;
874	MMcFadcHeader **HeaderFadc = NULL;
875
876	int numberBranches;
877
878	if (!Trigger_Loop)
879	numberBranches=ct_Number;
880	else
881	numberBranches=icontrigger;
882
883	if (Write_McTrig){
884
885	McTrig = new MMcTrig * [numberBranches];
886
887	for (int i=0;i<numberBranches;i++) {
888	McTrig[i] = new MMcTrig();
889	}
890
891	HeaderTrig = new MMcTrigHeader * [numberBranches];
892
893	for (int i=0;i<numberBranches;i++) {
894	HeaderTrig[i] = new MMcTrigHeader();
895	}
896	}
897
898	if (Write_McFADC){
899
900	HeaderFadc = new MMcFadcHeader * [numberBranches];
901
902	for (int i=0;i<numberBranches;i++) {
903	HeaderFadc[i] = new MMcFadcHeader();
904	}
905	}
906
907	MFadc **Fadc_CT;
908	Fadc_CT = new MFadc *[ct_Number];
909
910	for (int i=0; i<ct_Number;i++)
911	Fadc_CT[i] =
912	new MFadc(((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels(),
913	FADC_response_ampl,FADC_response_fwhm,
914	FADC_resp_ampl_out,FADC_resp_fwhm_out) ; //@< A instance of the Class MFadc
915
916
917	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
918	//
919	// Set the FADC pedestals for that run
920	// Some modifications
921	// mut be done to simulate a more realistic distribution of the pedestals.
922	// This simualtion is done int the SetPedestals methode inside the
923	// class MFadc
924	// Currentlly a input_pedestal array is declared with the pedestals.
925	// Thy can also be set randomly following a flat distribution.
926	//
927	/////////////////////////////////////////////////////////////////////
928
929	Float_t input_pedestals[ct_NPixels];
930
931	for(int i=0;i<ct_NPixels;i++)
932	input_pedestals[i]=get_FADC_pedestal();
933	for (int ict=0; ict<ct_Number;ict++){
934
935	Fadc_CT[ict]->SetPedestals(input_pedestals);
936	Fadc_CT[ict]->SetHigh2LowGain(FADC_high2low);
937
938	// Generate databse for the Fadc electronic noise
939	if (addElecNoise){
940	Fadc_CT[ict]->SetElecNoise(FADC_noise);
941	Fadc_CT[ict]->SetDigitalNoise(DIGITAL_noise);
942	}
943	}
944
945	// Prepare the raw data output
946
947	// Header Tree
948	MGeomCam *camdummy = new MGeomCam [ct_Number];
949	for (int ict=0; ict<ct_Number;ict++){
950	camdummy[ict]= ((MGeomCam*)(camgeom.UncheckedAt(ict)));
951	}
952	MRawRunHeader *RunHeader= new MRawRunHeader();
953	MMcRunHeader *McRunHeader = new MMcRunHeader();
954	MMcCorsikaRunHeader *McCorsikaRunHeader = new MMcCorsikaRunHeader();
955
956	MMcConfigRunHeader **McConfigRunHeader = NULL;
957	McConfigRunHeader = new MMcConfigRunHeader * [numberBranches];
958	for (int i=0;i<numberBranches;i++) {
959	McConfigRunHeader[i] = new MMcConfigRunHeader();
960	}
961
962	// Header branch
963
964	MRawEvtHeader **EvtHeader = NULL;
965
966	if (Write_RawEvt) {
967	EvtHeader = new MRawEvtHeader * [numberBranches];
968
969	for (int i=0;i<numberBranches;i++) {
970	EvtHeader[i] = new MRawEvtHeader();
971	}
972	}
973
974	// Data branch
975
976	MRawEvtData **EvtData = NULL; // Data branch
977
978	if (Write_RawEvt) {
979	EvtData = new MRawEvtData * [numberBranches];
980
981	for (int i=0;i<numberBranches;i++) {
982	EvtData[i] = new MRawEvtData();
983	EvtData[i]->Init(RunHeader); // We need thr RunHeader to read
984	// number of pixels
985	}
986	}
987
988	MMcEvt **McEvt = NULL;
989
990	if (Write_McEvt) {
991	if (ct_Number>1){
992	McEvt = new MMcEvt *[numberBranches];
993	for (int i=0;i<numberBranches;i++)
994	McEvt[i]=new MMcEvt();
995	}
996	else {
997	McEvt = new MMcEvt *[1];
998	McEvt[0]=new MMcEvt();
999	}
1000	}
1001	//
1002	// initalize a temporal ROOT file
1003	//
1004
1005	TFile outfile_temp ( rootname , "RECREATE" );
1006
1007	// create a Tree for the Header Event
1008	TTree HeaderTree("RunHeaders","Headers of Run");
1009
1010	// define branches of Header Tree
1011
1012	char help[4];
1013
1014	HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
1015	&RunHeader);
1016
1017	HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
1018	&McRunHeader);
1019
1020	HeaderTree.Branch("MMcCorsikaRunHeader","MMcCorsikaRunHeader",
1021	&McCorsikaRunHeader);
1022
1023	if(!Trigger_Loop && Write_McTrig && ct_Number==1){
1024
1025	HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader",
1026	&HeaderTrig[0]);
1027	}
1028	if (ct_Number==1){
1029	HeaderTree.Branch("MGeomCam",GeometryName[0],
1030	&camdummy[0]);
1031	HeaderTree.Branch("MMcConfigRunHeader","MMcConfigRunHeader",
1032	&McConfigRunHeader[0]);
1033	}
1034	else{
1035	char branchname[20];
1036	for (int ict=0; ict<ct_Number;ict++){
1037	sprintf(help,"%i",ict+1);
1038	strcpy (branchname, "MGeomCam;");
1039	strcat (branchname, & help[0]);
1040	strcat (branchname, ".");
1041	HeaderTree.Branch(branchname,GeometryName[ict],
1042	&camdummy[ict]);
1043	}
1044	for(int i=0;i<numberBranches;i++){
1045	sprintf(help,"%i",i+1);
1046	strcpy (branchname, "MMcConfigRunHeader;");
1047	strcat (branchname, & help[0]);
1048	strcat (branchname, ".");
1049	HeaderTree.Branch(branchname,"MMcConfigRunHeader",
1050	&McConfigRunHeader[i]);
1051	}
1052
1053	}
1054
1055	if ((Trigger_Loop \|\| ct_Number>1) && Write_McTrig){
1056	ibr=0;
1057	for(char branchname[20];ibr<numberBranches;ibr++){
1058
1059	sprintf(help,"%i",ibr+1);
1060	strcpy (branchname, "MMcTrigHeader;");
1061	strcat (branchname, & help[0]);
1062	strcat (branchname, ".");
1063	HeaderTree.Branch(branchname,"MMcTrigHeader",
1064	&HeaderTrig[ibr]);
1065	}
1066	}
1067
1068	if(!Trigger_Loop && Write_McFADC && ct_Number==1){
1069
1070	HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader",
1071	&HeaderFadc[0]);
1072	}
1073	if ((Trigger_Loop \|\| ct_Number>1) && Write_McFADC){
1074	ibr=0;
1075	for(char branchname[20];ibr<numberBranches;ibr++){
1076
1077	sprintf(help,"%i",ibr+1);
1078	strcpy (branchname, "MMcFadcHeader;");
1079	strcat (branchname, & help[0]);
1080	strcat (branchname, ".");
1081	HeaderTree.Branch(branchname,"MMcFadcHeader",
1082	&HeaderFadc[ibr]);
1083	}
1084	}
1085
1086	// Fill branches for MRawRunHeader
1087
1088	RunHeader->SetMagicNumber(kMagicNumber);
1089	RunHeader->SetFormatVersion(4);
1090	RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
1091	RunHeader->SetRunType(256);
1092	RunHeader->SetRunNumber(0);
1093	RunHeader->SetNumSamples(0,FADC_SLICES);
1094	RunHeader->SetNumCrates(1);
1095	RunHeader->SetNumPixInCrate(ct_NPixels);
1096
1097	// Fill branches for MMcTrigHeader
1098
1099	if(!Trigger_Loop && Write_McTrig && ct_Number==1){
1100
1101	HeaderTrig[0]->SetTopology((Short_t) Trigger_topology[0]);
1102	HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity[0]);
1103	HeaderTrig[0]->SetThreshold(qThreshold[0]);
1104	HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
1105	HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
1106	HeaderTrig[0]->SetOverlap(Trigger_overlaping_time);
1107	HeaderTrig[0]->SetGate(Trigger_gate_length);
1108	HeaderTrig[0]->SetElecNoise(Trigger_noise);
1109	}
1110	if(!Trigger_Loop && Write_McTrig && ct_Number>1){
1111	for(int i=0;i<ct_Number;i++){
1112	HeaderTrig[i]->SetTopology((Short_t) Trigger_topology[i]);
1113	HeaderTrig[i]->SetMultiplicity((Short_t) Trigger_multiplicity[i]);
1114	HeaderTrig[i]->SetThreshold(qThreshold[i]);
1115	HeaderTrig[i]->SetAmplitud(Trigger_response_ampl);
1116	HeaderTrig[i]->SetFwhm(Trigger_response_fwhm);
1117	HeaderTrig[i]->SetOverlap(Trigger_overlaping_time);
1118	HeaderTrig[i]->SetGate(Trigger_gate_length);
1119	HeaderTrig[i]->SetElecNoise(Trigger_noise);
1120	}
1121	}
1122	if(Trigger_Loop && Write_McTrig){
1123
1124	int iconcount;
1125	for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++,fthrescount+=Trigger_loop_sthres){
1126	for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
1127	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
1128	HeaderTrig[iconcount]->SetTopology((Short_t) isorttopo[itopocount+Trigger_loop_ltop]);
1129	HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
1130	for(int i=0;i<ct_NPixels;i++){
1131	fpixelthres[i]=
1132	((Float_t)(fthrescount)>=qThreshold[0][i])?
1133	(Float_t)(fthrescount):qThreshold[0][i];
1134	}
1135	HeaderTrig[iconcount]->SetThreshold( fpixelthres);
1136	HeaderTrig[iconcount]->SetAmplitud(Trigger_response_ampl);
1137	HeaderTrig[iconcount]->SetFwhm(Trigger_response_fwhm);
1138	HeaderTrig[iconcount]->SetOverlap(Trigger_overlaping_time);
1139	HeaderTrig[iconcount]->SetGate(Trigger_gate_length);
1140	HeaderTrig[iconcount]->SetElecNoise(Trigger_noise);
1141	iconcount++;
1142	}
1143	}
1144	}
1145	}
1146
1147	// Fill branches for MMcFadcHeader
1148
1149	for(int i=0;i<ct_NPixels;i++){
1150	fadc_elecnoise[i]=FADC_noise;
1151	fadc_diginoise[i]=DIGITAL_noise;
1152	}
1153
1154	Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
1155
1156	if(!Trigger_Loop && Write_McFADC && ct_Number==1){
1157
1158	HeaderFadc[0]->SetShape(0.0);
1159	HeaderFadc[0]->SetAmplitud(FADC_response_ampl*2.35/
1160	sqrt(2M_PIFADC_response_fwhm
1161	*FADC_response_fwhm),
1162	FADC_resp_ampl_out*2.35/
1163	sqrt(2M_PIFADC_resp_fwhm_out
1164	*FADC_resp_fwhm_out));
1165	HeaderFadc[0]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
1166	HeaderFadc[0]->SetLow2High(FADC_high2low);
1167	HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],
1168	((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
1169	HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
1170	((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
1171	}
1172	if(!Trigger_Loop && Write_McFADC && ct_Number>1){
1173	for(int i=0;i<ct_Number;i++){
1174	Fadc_CT[i]->GetPedestals(&fadc_pedestals[0]);
1175	HeaderFadc[i]->SetShape(0.0);
1176	HeaderFadc[i]->SetAmplitud(FADC_response_ampl*2.35/
1177	sqrt(2M_PIFADC_response_fwhm
1178	*FADC_response_fwhm),
1179	FADC_resp_ampl_out*2.35/
1180	sqrt(2M_PIFADC_resp_fwhm_out
1181	*FADC_resp_fwhm_out));
1182	HeaderFadc[i]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
1183	HeaderFadc[i]->SetLow2High(FADC_high2low);
1184	HeaderFadc[i]->SetPedestal(&fadc_pedestals[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
1185	HeaderFadc[i]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
1186	((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
1187	}
1188	}
1189	if(Trigger_Loop && Write_McFADC){
1190	int iconcount;
1191	for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1192	for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
1193	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
1194	Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
1195	HeaderFadc[iconcount]->SetShape(0.0);
1196	HeaderFadc[iconcount]->SetAmplitud(FADC_response_ampl*2.35/
1197	sqrt(2M_PIFADC_response_fwhm
1198	*FADC_response_fwhm),
1199	FADC_resp_ampl_out*2.35/
1200	sqrt(2M_PIFADC_resp_fwhm_out
1201	*FADC_resp_fwhm_out));
1202	HeaderFadc[iconcount]->SetFwhm(FADC_response_fwhm,
1203	FADC_resp_fwhm_out);
1204	HeaderFadc[iconcount]->SetLow2High(FADC_high2low);
1205	HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0], ct_NPixels);
1206	HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],
1207	&fadc_diginoise[0],ct_NPixels);
1208	iconcount++;
1209	}
1210	}
1211	}
1212	}
1213
1214	// UnSet flag in pixel 0 (not used for trigger) that indicates if secure pixel
1215	// is active once it has been stored
1216
1217	for (int ict=0; ict<ct_Number;ict++){
1218	qThreshold[ict][0]=999999.99;
1219	}
1220
1221
1222	// create a Tree for the Event data stream
1223	TTree EvtTree("Events","Normal Triggered Events");
1224
1225	if (Write_McEvt && ct_Number==1){
1226
1227	EvtTree.Branch("MMcEvt","MMcEvt",
1228	&McEvt[0]);
1229	}
1230
1231	if (Write_McEvt && ct_Number!=1){
1232	ibr=0;
1233	for(char branchname[10];ibr<numberBranches;ibr++){
1234
1235	sprintf(help,"%i",ibr+1);
1236	strcpy (branchname, "MMcEvt;");
1237	strcat (branchname, & help[0]);
1238	strcat (branchname, ".");
1239	EvtTree.Branch(branchname,"MMcEvt",
1240	&McEvt[ibr]);
1241	}
1242	}
1243
1244	if(!Trigger_Loop && ct_Number==1){
1245
1246	if (Write_RawEvt){
1247	EvtTree.Branch("MRawEvtHeader","MRawEvtHeader",
1248	&EvtHeader[0]);
1249	EvtTree.Branch("MRawEvtData","MRawEvtData",
1250	&EvtData[0]);
1251	}
1252	if (Write_McTrig){
1253	EvtTree.Branch("MMcTrig","MMcTrig",
1254	&McTrig[0]);
1255	}
1256	}
1257	else{
1258
1259	if (Write_McTrig){
1260	ibr=0;
1261	for(char branchname[10];ibr<numberBranches;ibr++){
1262
1263	sprintf(help,"%i",ibr+1);
1264	strcpy (branchname, "MMcTrig;");
1265	strcat (branchname, & help[0]);
1266	strcat (branchname, ".");
1267	EvtTree.Branch(branchname,"MMcTrig",
1268	&McTrig[ibr]);
1269	}
1270	}
1271	}
1272
1273	if ((Trigger_Loop \|\| ct_Number>1) && Write_RawEvt){
1274	ibr=0;
1275	for(char branchname[15];ibr<numberBranches;ibr++){
1276
1277	sprintf(help,"%i",ibr+1);
1278	strcpy (branchname, "MRawEvtHeader;");
1279	strcat (branchname, & help[0]);
1280	strcat (branchname, ".");
1281	EvtTree.Branch(branchname,"MRawEvtHeader",
1282	&EvtHeader[ibr]);
1283	}
1284	ibr=0;
1285	for(char branchname[15];ibr<numberBranches;ibr++){
1286
1287	sprintf(help,"%i",ibr+1);
1288	strcpy (branchname, "MRawEvtData;");
1289	strcat (branchname, & help[0]);
1290	strcat (branchname, ".");
1291	EvtTree.Branch(branchname,"MRawEvtData",
1292	&EvtData[ibr]);
1293	}
1294	}
1295
1296
1297	TApplication theAppTrigger("App", &argc, argv);
1298
1299	if(FADC_Scan){
1300	if (gROOT->IsBatch()) {
1301	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1302	// return 1;
1303	}
1304	}
1305	if(FADC_Scan){
1306	//TApplication theAppFadc("App", &argc, argv);
1307
1308	if (gROOT->IsBatch()) {
1309	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1310	// return 1;
1311	}
1312	}
1313
1314
1315	// prepare the NSB simulation
1316
1317	// Instance of the Mlons class
1318	MLons lons(Trigger_response_ampl, Trigger_response_fwhm,
1319	FADC_response_ampl,FADC_response_fwhm);
1320
1321	lons.SetPath(nsbpathname);
1322
1323	// Instance of the Mlons class
1324	MLons lons_outer(Trigger_response_ampl, Trigger_response_fwhm,
1325	FADC_resp_ampl_out,FADC_resp_fwhm_out);
1326
1327	lons_outer.SetPath(nsbpath_outer);
1328
1329	if( simulateNSB ){
1330
1331	//
1332	// Calculate the non-diffuse NSB photoelectron rates
1333	//
1334
1335	// FIXME --- star NSB different for each camera?
1336	log(SIGNATURE,"Produce NSB rates from Star Field");
1337
1338	k = produce_nsbrates( starfieldname,
1339	((MGeomCam*)(camgeom.UncheckedAt(0))),
1340	nsbrate_phepns,
1341	0);
1342
1343	if (k != 0){
1344	cout << "Error when reading starfield... \nExiting.\n";
1345	exit(1);
1346	}
1347
1348	// calculate diffuse rate correcting for the pixel size and telescope
1349
1350	float factorNSB_ct;
1351	for(int ict=0;ict<ct_Number;ict++){
1352
1353	// First we set the factor due to the mirror size respect to the MAGIC mirror.
1354
1355	switch(GeometryCamera[ict]){
1356	case 1:
1357	case 2:
1358	case 3:
1359	case 4:
1360	factorNSB_ct=1.0;
1361	break;
1362	case 5:
1363	case 6:
1364	case 7:
1365	factorNSB_ct=1000.0/239.0;
1366	break;
1367	case 8:
1368	case 9:
1369	default:
1370	factorNSB_ct=1.0;
1371	break;
1372	}
1373
1374	factorNSB_ct=factorNSB_ct/
1375	((((MGeomCam)(camgeom.UncheckedAt(ict)))).GetCameraDist()1000
1376	(((MGeomCam)(camgeom.UncheckedAt(ict)))).GetCameraDist()1000
1377	PIXEL_SIZE*PIXEL_SIZE);
1378
1379	for(UInt_t ui=0;
1380	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); ui++){
1381	const Float_t size=
1382	(((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD();
1383	diffnsb_phepns[ict][ui] = (Int_t(meanNSBfactorNSB_ct100sizesize+0.5))/(100.0)*factorqe_NSB[ict];
1384	}
1385	}
1386
1387	// calculate nsb rate including diffuse and starlight
1388	// we also include the extinction effect
1389	for(int ict=0;ict<ct_Number;ict++){
1390	for(j=0;j<iNUMWAVEBANDS;j++){
1391	// calculate the effect of the atmospheric extinction
1392
1393	zenfactor = pow(10., -0.4 * ext[j] );
1394
1395	for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1396	nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
1397	nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
1398	}
1399	}
1400	}
1401	}
1402
1403
1404	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1405	//
1406	// Now 500 empty event with the conditin in which the camera is run
1407	// are simulated. IN this way one gets an estimation of the
1408	// sigma for the pedestal of each FADC channel.
1409	// This sigma computtaion is done assuming any noise taht effects
1410	// the FADC but there is no rotation of the Star Field (otherwise it
1411	// should be done for each event)
1412	//
1413	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1414
1415	for(int ict=0;ict<ct_Number;ict++){
1416	for(UInt_t ui=0;
1417	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1418	ui++){
1419	fadc_sigma[ui]=0.0;
1420	fadc_sigma_low[ui]=0.0;
1421	}
1422	for(int ie=0;ie<500;ie++){
1423	Fadc_CT[ict]->Reset() ;
1424	if (addElecNoise){
1425	Fadc_CT[ict]->ElecNoise(FADC_noise) ;
1426	}
1427	if(simulateNSB){
1428	for(UInt_t ui=0;
1429	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1430	ui++){
1431	if(nsb_phepns[ict][ui]>0.0){
1432	if((((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD()>(((MGeomCam)(camgeom.UncheckedAt(ict))))[0].GetD()){
1433	k=lons_outer.GetResponse(nsb_phepns[ict][ui],0.01,
1434	& nsb_trigresp[0],
1435	& nsb_fadcresp[0]);
1436	}
1437	else{
1438	k=lons.GetResponse(nsb_phepns[ict][ui],0.01,
1439	& nsb_trigresp[0],& nsb_fadcresp[0]);
1440	}
1441	if(k==0){
1442	cout << "Exiting.\n";
1443	exit(1);
1444	}
1445	Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
1446	}
1447	}
1448	}
1449	Fadc_CT[ict]->Pedestals();
1450	for(UInt_t ui=0;
1451	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1452	ui++){
1453	fadc_sigma[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,1))/500.0;
1454	fadc_sigma_low[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,0))/500.0;
1455	}
1456	}
1457	HeaderFadc[ict]->SetPedestalSigma(&fadc_sigma_low[0],&fadc_sigma[0],
1458	((MGeomCam*)(camgeom.UncheckedAt(ict)))
1459	->GetNumPixels());
1460	}
1461
1462	//
1463	// Read the reflector file with the Cherenkov data
1464	//
1465
1466	// select input file
1467
1468	if ( Data_From_STDIN ) {
1469
1470	inputfile[0] = stdin;
1471
1472	}
1473	else{
1474
1475	for(int ict=0;ict<ct_Number;ict++){
1476
1477	log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
1478	inputfile[ict] = fopen( inname_CT[ict], "r" );
1479
1480	if ( inputfile[ict] == NULL )
1481	error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
1482	}
1483
1484	}
1485
1486	// get signature, and check it
1487
1488	for(int ict=0;ict<ct_Number;ict++){
1489	if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE){
1490	exit(1);
1491	}
1492	}
1493
1494	// open data file
1495
1496	log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
1497	datafile.open( datname );
1498
1499	if ( datafile.bad() )
1500	error( SIGNATURE, "Cannot open data file: %s\n", datname );
1501
1502	// initializes flag
1503
1504	strcpy( flag, " \0" );
1505	strcpy( flag_new, " \0" );
1506
1507	// allocate space for PMTs numbers of pixels
1508
1509	fnpix = new float [ct_NPixels];
1510
1511
1512	//!@' @#### Main loop.
1513	//@'
1514
1515	// get flag
1516
1517	for(int ict=0;ict<ct_Number;ict++)
1518	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1519
1520
1521	// loop over the file
1522
1523	still_in_loop = TRUE;
1524
1525	// FIXME --- check if not eof for all input reflector files
1526
1527	while (
1528	((! Data_From_STDIN) && ( !feof(inputfile[0]) ))
1529	\|\|
1530	(Data_From_STDIN && still_in_loop)
1531	) {
1532
1533	// reading .rfl files
1534	if(!isA( flag, FLAG_START_OF_RUN )){
1535
1536	// We break the main loop
1537	cout<<"Warning: Expected start of run flag, but found:"<<flag<<endl;
1538	cout<<" We break the main loop"<<endl;
1539	break;
1540	}
1541	else { // found start of run
1542
1543	for(int ict=0;ict<ct_Number;ict++) {
1544
1545	fread( flag_new, 4, 1, inputfile[ict] );
1546
1547	if(!isA( flag_new, FLAG_START_OF_HEADER)){
1548
1549	// We break the main loop
1550	cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
1551	cout<<" We break the main loop"<<endl;
1552	break;
1553	}
1554
1555	fread((char)&mcrunh,(SIZE_OF_MCRUNHEADER)sizeof(float), 1, inputfile[ict] );
1556	}
1557
1558	nshow=0;
1559
1560	for(int ict=0;ict<ct_Number;ict++){
1561	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1562	}
1563	while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
1564	for(int ict=0;ict<ct_Number;ict++) {
1565	fread( flag_new, 4, 1, inputfile[ict] );
1566	if(!isA( flag_new, FLAG_EVENT_HEADER)){
1567
1568	// We break while events loop
1569	cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<" "<<ict<<endl;
1570	cout<<" We break the while events loop"<<endl;
1571	break;
1572	}
1573	}
1574
1575
1576	//
1577	// Clear Trigger and Fadc
1578	//
1579	for(int ict=0;ict<ct_Number;ict++) {
1580	Trigger_CT[ict]->Reset() ;
1581	Trigger_CT[ict]->ClearFirst();
1582	Trigger_CT[ict]->ClearZero();
1583	Fadc_CT[ict]->Reset() ;
1584	}
1585
1586	++nshow;
1587	if((nshow+ntshow)%100 == 1)
1588	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
1589
1590	// get MCEventHeader
1591	if (reflector_file_version<6){
1592	for(int ict=0;ict<ct_Number;ict++)
1593	fread( (char*)&mcevth[ict], mcevth[ict].mysize(),
1594	1, inputfile[ict] );
1595	}
1596	else{
1597	for(int ict=0;ict<ct_Number;ict++)
1598	fread( (char*)&mcevth_2[ict], mcevth_2[ict].mysize(),
1599	1, inputfile[ict] );
1600	}
1601
1602	// calculate impact parameter (shortest distance betwee the original
1603	// trajectory of the primary (assumed shower-axis) and the
1604	// direction where the telescope points to
1605	//
1606	// we use the following equation, given that the shower core position
1607	// is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
1608	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
1609	// and (l2,m2,n2)
1610	//
1611	// \| \|
1612	// \| x1-x2 y1-y2 z1-z2 \|
1613	// \| \|
1614	// + \| l1 m1 n1 \|
1615	// - \| \|
1616	// \| l2 m2 n2 \|
1617	// \| \|
1618	// dist = ------------------------------------ ( > 0 )
1619	// [ \|l1 m1\|2 \|m1 n1\|2 \|n1 l1\|2 ]1/2
1620	// [ \| \| + \| \| + \| \| ]
1621	// [ \|l2 m2\| \|m2 n2\| \|n2 l2\| ]
1622	//
1623	// playing a little bit, we get this reduced for in our case:
1624	//
1625	//
1626	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
1627	// [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
1628	// 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
1629
1630	// read the direction of the incoming shower
1631	// It is done only for one telescope since it is suposed
1632	// to be the same shower for all of them
1633	if (reflector_file_version<6){
1634	thetashw = mcevth[0].get_theta();
1635	phishw = mcevth[0].get_phi();
1636	}
1637	else{
1638	thetashw = mcevth_2[0].get_theta();
1639	phishw = mcevth_2[0].get_phi();
1640	}
1641
1642	// calculate vector for shower
1643
1644	l1 = sin(thetashw)*cos(phishw);
1645	m1 = sin(thetashw)*sin(phishw);
1646	n1 = cos(thetashw);
1647
1648	if (reflector_file_version<6){
1649
1650	mcevth[0].get_core(&coreX, &coreY);
1651
1652	for(int ict=0;ict<ct_Number;ict++){
1653	// read the deviation of the telescope with respect to the shower
1654	mcevth[ict].get_deviations ( &thetaCT[ict], &phiCT[ict] );
1655
1656	if ( (thetaCT[ict] == 0.) && (phiCT[ict] == 0.) ) {
1657
1658	// CT was looking to the source (both lines are parallel)
1659	// therefore, we calculate the impact parameter as the distance
1660	// between the CT axis and the core position
1661
1662	impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1663	thetaCT[ict] += thetashw;
1664	phiCT[ict] += phishw;
1665	} else {
1666
1667	// the shower comes off-axis
1668
1669	// obtain with this the final direction of the CT
1670
1671	thetaCT[ict] += thetashw;
1672	phiCT[ict] += phishw;
1673
1674	// calculate vector for telescope
1675
1676	l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
1677	m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
1678	n2 = cos(thetaCT[ict]);
1679
1680	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1681	den = (SQR(l1m2 - l2m1) +
1682	SQR(m1n2 - m2n1) +
1683	SQR(n1l2 - n2l1));
1684	den = sqrt(den);
1685	impactD[ict] = fabs(num)/den;
1686	}
1687	}
1688	}
1689	else{
1690	mcevth_2[0].get_core(&coreX, &coreY);
1691	for(int ict=0;ict<ct_Number;ict++){
1692	thetaCT[ict]=mcevth_2[ict].get_theta_CT();
1693	phiCT[ict]=mcevth_2[ict].get_phi_CT();
1694	if ( (thetaCT[ict] == thetashw) && (phiCT[ict] == phishw) ) {
1695
1696	// CT was looking to the source (both lines are parallel)
1697	// therefore, we calculate the impact parameter as the distance
1698	// between the CT axis and the core position
1699
1700	impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1701
1702	} else {
1703
1704	// the shower comes off-axis
1705
1706	// calculate vector for telescope
1707
1708	l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
1709	m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
1710	n2 = cos(thetaCT[ict]);
1711
1712	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1713	den = (SQR(l1m2 - l2m1) +
1714	SQR(m1n2 - m2n1) +
1715	SQR(n1l2 - n2l1));
1716	den = sqrt(den);
1717
1718	impactD[ict] = fabs(num)/den;
1719
1720	}
1721	}
1722	}
1723
1724	// energy cut
1725
1726	if ( Select_Energy ) {
1727	if (reflector_file_version<6)
1728	if (( mcevth[0].get_energy() < Select_Energy_le ) \|\|
1729	( mcevth[0].get_energy() > Select_Energy_ue )) {
1730	log(SIGNATURE, "select_energy: shower rejected.\n");
1731	continue;
1732	}
1733	else
1734	if (( mcevth_2[0].get_energy() < Select_Energy_le ) \|\|
1735	( mcevth_2[0].get_energy() > Select_Energy_ue )) {
1736	log(SIGNATURE, "select_energy: shower rejected.\n");
1737	continue;
1738	}
1739	}
1740
1741	// Read first and last time and put inumphe_CT[0] to 0
1742	float tm=0, tM=9E+50;
1743	for(int ict=0;ict<ct_Number;ict++){
1744
1745	if (reflector_file_version<6)
1746	mcevth[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1747	else
1748	mcevth_2[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1749
1750	tm=(tm<arrtmin_ns)?tm:arrtmin_ns;
1751	tM=(tM<arrtmax_ns)?tM:arrtmax_ns;
1752
1753	}
1754	arrtmin_ns=tm;arrtmax_ns=tM;
1755
1756	inumphe=0;
1757
1758	for(int ict=0;ict<ct_Number;ict++){
1759
1760	inumphe_CT[ict]=0;
1761
1762	// read the photons and produce the photoelectrons
1763
1764	k = produce_phes( inputfile[ict],
1765	((MGeomCam*)(camgeom.UncheckedAt(ict))),
1766	WAVEBANDBOUND1,
1767	WAVEBANDBOUND6,
1768	Trigger_CT[ict], // will be changed by the function!
1769	Fadc_CT[ict], // will be changed by the function!
1770	&inumphe_CT[ict], // important for later: the size of photoe[]
1771	fnpix, // will be changed by the function!
1772	&ncph[ict], // will be changed by the function!
1773	&arrtmin_ns, // will be changed by the function!
1774	&arrtmax_ns, // will be changed by the function!
1775	ict);
1776
1777	inumphe=(inumphe<inumphe_CT[ict])?inumphe_CT[ict]:inumphe;
1778
1779	if( k != 0 ){ // non-zero returnvalue means error
1780	cout << "Exiting.\n";
1781	exit(1);
1782	}
1783	}
1784
1785	// NSB simulation
1786
1787	if(simulateNSB && nphe2NSB<=inumphe){
1788
1789	if(Starfield_rotate){
1790
1791	// Introduction rho angle
1792
1793	zenith = thetashw;
1794	azimutal = phishw;
1795	C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
1796	C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
1797	C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
1798	rho = acos( C1/C2 );
1799
1800	if ( sin(azimutal) < 0)
1801	{
1802	rho = 2 * 3.14159 - rho;
1803	}
1804	else
1805	{
1806	rho = rho;
1807	}
1808	rho = rho*180/3.14159;
1809
1810	// Rottion of the NSB
1811	// FIXME --- We should rotate for all cameras. Is it always the same rho?
1812	for(int ict=0;ict<ct_Number;ict++){
1813	k = size_rotated(
1814	&nsb_phepns_rotated[ict][0],
1815	nsb_phepns[ict],
1816	rho);
1817	}
1818	}
1819
1820	// Fill trigger and fadc response in the trigger class from the database
1821	for(int ict=0;ict<ct_Number;ict++){
1822
1823	for(UInt_t ui=0;ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1824	if(nsb_phepns_rotated[ict][ui]>0.0){
1825	if((((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD()>(((MGeomCam)(camgeom.UncheckedAt(ict))))[0].GetD()){
1826	k=lons_outer.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
1827	& nsb_trigresp[0],
1828	& nsb_fadcresp[0]);
1829	}
1830	else{
1831	k=lons.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
1832	& nsb_trigresp[0],& nsb_fadcresp[0]);
1833	}
1834	if(k==0){
1835	cout << "Exiting.\n";
1836	exit(1);
1837	}
1838	Trigger_CT[ict]->AddNSB(ui,nsb_trigresp);
1839	Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
1840	}
1841	}
1842	}
1843
1844	}// end if(simulateNSB && nphe2NSB<=inumphe_CT[0]) ...
1845
1846	for(int ict=0;ict<ct_Number;ict++){
1847	inumphensb[ict]=0;
1848	for (UInt_t ui=0;ui<
1849	((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1850	inumphensb[ict]+=nsb_phepns[ict][ui]*TOTAL_TRIGGER_TIME;
1851	}
1852	ntcph[ict]+=ncph[ict];
1853	if ((nshow+ntshow)%100 == 1){
1854	log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1855	ncph[ict], ntcph[ict]);
1856
1857	cout << "Total number of phes in CT "<<ict<<": "
1858	<< inumphe_CT[ict]<<" (+ ";
1859	cout<<inumphensb[ict]<<" mean expected number from NSB)"<<endl;
1860	}
1861	}
1862
1863	// skip it ?
1864
1865	int i;
1866	for (i=0; i<nSkip; ++i ) {
1867	if (Skip[i] == (nshow+ntshow)) {
1868	i = -1;
1869	break;
1870	}
1871	}
1872
1873	// if after the previous loop, the exit value of i is -1
1874	// then the shower number is in the list of showers to be
1875	// skipped
1876
1877	if (i == -1) {
1878	log(SIGNATURE, "\t\tskipped!\n");
1879	continue;
1880	}
1881
1882	//++++++++++++++++++++++++++++++++++++++++++++++++++
1883	// at this point we have a camera full of
1884	// ph.e.s
1885	// we should first apply the trigger condition,
1886	// and if there's trigger, then clean the image,
1887	// calculate the islands statistics and the
1888	// other parameters of the image (Hillas' parameters
1889	// and so on).
1890	//--------------------------------------------------
1891
1892	// TRIGGER HERE
1893
1894	// We should simulate the AC coupling behaviour:
1895	// For the FADC it is only done for the NSB while producinr
1896	// the StarResponse database.
1897	// For the trigger is done in the Trigger.Diskriminate(), which
1898	// is called later (it should be seperated to speed up the program
1899	//
1900
1901	// now the noise of the electronic
1902	// (preamps, optical transmission,..) is introduced.
1903	// This is done inside the class MTrigger by the method ElecNoise.
1904	//
1905
1906	for(int ict=0;ict<ct_Number;ict++) {
1907	if (addElecNoise && nphe2NSB<=inumphe){
1908
1909	Trigger_CT[ict]->ElecNoise(Trigger_noise) ;
1910
1911	Fadc_CT[ict]->ElecNoise(FADC_noise) ;
1912	}
1913	}
1914
1915	// now a shift in the fadc signal due to the pedestlas is
1916	// intorduced
1917	// This is done inside the class MFadc by the method Pedestals
1918	for(int ict=0;ict<ct_Number;ict++) {
1919	Fadc_CT[ict]->Pedestals();
1920	}
1921
1922	// We study several trigger conditons
1923	if(Trigger_Loop){
1924
1925	// Set to zero the flag to know if some conditon has triggered
1926	btrigger=0;
1927	flagstoring = 0;
1928
1929	// Loop over trigger threshold
1930	int iconcount;
1931	for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1932	for (int i=0;i<ct_NPixels;i++){
1933	fpixelthres[i]=
1934	((Float_t)(fthrescount)>=qThreshold[0][i])?
1935	(Float_t)(fthrescount):qThreshold[0][i];
1936
1937	// Rise the discrimnator threshold to avoid huge rates
1938
1939	if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
1940	for(int ii=0;ii<ct_NPixels;ii++)
1941	if( nsb_phepns_rotated[0][ii]>riseDiskThres)
1942	fpixelthres[ii]=secureDiskThres;
1943	}
1944	Trigger_CT[0]->SetThreshold(fpixelthres);
1945
1946	Trigger_CT[0]->Diskriminate();
1947
1948
1949	//
1950	// look if in all the signals in the trigger signal branch
1951	// is a possible Trigger. Therefore we habe to diskriminate all
1952	// the simulated analog signals (Method Diskriminate in class
1953	// MTrigger). We look simultanously for the moments at which
1954	// there are more than TRIGGER_MULTI pixels above the
1955	// CHANNEL_THRESHOLD.
1956	//
1957
1958	// Set trigger flags to zero
1959	Lev0=0;
1960	Lev1=0;
1961
1962	// loop over multiplicity of trigger configuration
1963	for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1964	Trigger_CT[0]->SetMultiplicity(imulticount);
1965	Trigger_CT[0]->ClearZero();
1966
1967	Lev0=(Short_t) Trigger_CT[0]->ZeroLevel();
1968	if (Lev0>0 \|\| Write_All_Images \|\| btrigger){
1969
1970	// loop over topologies
1971	for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1972	Lev1=0;
1973
1974	if(itopocount==0 && imulticount>7) continue;
1975	//COBB if(itopocount==2 && imulticount<3) continue;
1976	// It only makes to look for a different topology
1977	// if there are 3 or more N pixels.
1978	if(imulticount<3)
1979	Trigger_CT[0]->SetTopology(1);
1980	else
1981	{
1982	// We should be careful that topologies are sort from
1983	// the less to the more restrictive one.
1984	Trigger_CT[0]->SetTopology(isorttopo[itopocount]);
1985	}
1986	Trigger_CT[0]->ClearFirst();
1987
1988	//
1989	// Start the First Level Trigger simulation
1990	//
1991	if(Lev0!=0)
1992	Lev1=Trigger_CT[0]->FirstLevel();
1993	if(Lev1>0) {
1994	btrigger= 1;
1995	ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++;
1996	}
1997
1998	Lev0=1;
1999	Int_t NumImages = Lev1;
2000	if(Lev1==0 && (Write_All_Images \|\| btrigger)){
2001	btrigger= 1;
2002	NumImages=1;
2003	Lev0=0;
2004	}
2005
2006	for (Int_t ii=0;ii<NumImages;ii++){
2007	if (Write_McTrig){
2008	McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
2009	McTrig[iconcount]->SetTime(Trigger_CT[0]->GetFirstLevelTime(ii),ii+1);
2010	Trigger_CT[0]->GetMapDiskriminator(trigger_map);
2011	McTrig[iconcount]->SetMapPixels(trigger_map,ii);
2012	}
2013	//
2014	// fill inside the class fadc the member output
2015	//
2016
2017	Fadc_CT[0]->TriggeredFadc(Trigger_CT[0]->GetFirstLevelTime(ii));
2018
2019	if( Write_RawEvt ){
2020	//
2021	// Fill the header of this event
2022	//
2023
2024	EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
2025	// fill pixel information
2026	if (Lev1 \|\| Write_All_Images){
2027	if (addElecNoise) Fadc_CT[0]->DigitalNoise();
2028	for(UInt_t i=0;
2029	i<((MGeomCam*)(camgeom.UncheckedAt(0)))
2030	->GetNumPixels();i++){
2031	if(!Fadc_CT[0]->IsPixelUsed(i)) continue;
2032	for (j=0;j<FADC_SLICES;j++){
2033	fadcValues->AddAt(Fadc_CT[0]->GetFadcSignal(i,j),j);
2034	fadcValuesLow->AddAt(Fadc_CT[0]->GetFadcLowGainSignal(i,j),j);
2035	}
2036	EvtData[iconcount]->AddPixel(i,fadcValues,0);
2037	EvtData[iconcount]->AddPixel(i,fadcValuesLow,kTRUE);
2038	}
2039	}
2040	}
2041	}
2042	//
2043	// Increase counter of analised trigger conditions
2044	//
2045	iconcount++;
2046	}
2047	}
2048	else{
2049	break;
2050	}
2051	}
2052	if (!btrigger) break;
2053	}
2054	if (btrigger){
2055
2056	//
2057	// fill the MMcEvt with all information
2058	//
2059
2060	if(!flagstoring)
2061	nstoredevents++;
2062	flagstoring = 1;
2063
2064	if (Write_McEvt) {
2065	Float_t ftime, ltime;
2066	if (reflector_file_version<6){
2067	mcevth[0].get_times(&ftime, &ltime);
2068	McEvt[0]->Fill( 0,
2069	(UShort_t) mcevth[0].get_primary() ,
2070	mcevth[0].get_energy(),
2071	-1.0,
2072	-1.0,
2073	-1.0,
2074	mcevth[0].get_theta(),
2075	mcevth[0].get_phi(),
2076	mcevth[0].get_core(),
2077	mcevth[0].get_coreX(),
2078	mcevth[0].get_coreY(),
2079	impactD[0],
2080	phiCT[0],
2081	thetaCT[0],
2082	ftime,
2083	ltime,
2084	0,
2085	0,
2086	0,
2087	0,
2088	0,
2089	0,
2090	0,
2091	(UInt_t)mcevth[0].get_CORSIKA(),
2092	(UInt_t)mcevth[0].get_AtmAbs(),
2093	(UInt_t)(mcevth[0].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()),
2094	(UInt_t) ncph[0],
2095	(UInt_t) inumphe_CT[0],
2096	(UInt_t) inumphensb[0]+inumphe_CT[0],
2097	-1.0,
2098	-1.0,
2099	-1.0);
2100	}
2101	else{
2102	Float_t Nmax, t0, tmax, a, b, c, chi2;
2103	mcevth_2[0].get_times(&ftime, &ltime);
2104	chi2=mcevth_2[0].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2105	McEvt[0]->Fill((UInt_t) mcevth_2[0].get_evt_number(),
2106	(UShort_t) mcevth_2[0].get_primary() ,
2107	mcevth_2[0].get_energy(),
2108	mcevth_2[0].get_thick0(),
2109	mcevth_2[0].get_first_target(),
2110	mcevth_2[0].get_z_first_int(),
2111	mcevth_2[0].get_theta(),
2112	mcevth_2[0].get_phi(),
2113	mcevth_2[0].get_core(),
2114	mcevth_2[0].get_coreX(),
2115	mcevth_2[0].get_coreY(),
2116	impactD[0],
2117	mcevth_2[0].get_phi_CT(),
2118	mcevth_2[0].get_theta_CT(),
2119	ftime,
2120	ltime,
2121	Nmax,
2122	t0,
2123	tmax,
2124	a,
2125	b,
2126	c,
2127	chi2,
2128	(UInt_t)mcevth_2[0].get_CORSIKA(),
2129	(UInt_t)mcevth_2[0].get_AtmAbs(),
2130	(UInt_t)(mcevth_2[0].get_MirrAbs()+mcevth_2[0].get_OutOfMirr()+mcevth_2[0].get_BlackSpot()),
2131	(UInt_t) ncph[0],
2132	(UInt_t) inumphe_CT[0],
2133	(UInt_t) inumphensb[0]+inumphe_CT[0],
2134	mcevth_2[0].get_ElecFraction(),
2135	mcevth_2[0].get_MuonFraction(),
2136	mcevth_2[0].get_OtherFraction());
2137	}
2138	}
2139	// Fill the Tree with the current leaves of each branch
2140	i=EvtTree.Fill() ;
2141
2142	// Clear the branches
2143	if(Write_McTrig){
2144	for(int i=0;i<numberBranches;i++){
2145	McTrig[i]->Clear() ;
2146	}
2147	}
2148	if( Write_RawEvt ){
2149	for(int i=0;i<numberBranches;i++){
2150	EvtHeader[i]->Clear() ;
2151	EvtData[i]->DeletePixels();
2152	}
2153	}
2154	if (Write_McEvt)
2155	McEvt[0]->Clear() ;
2156	}
2157	}
2158	// We study a single trigger condition
2159	else {
2160
2161	// Set to zero the flag to know if some conditon has triggered
2162	btrigger=0;
2163	flagstoring = 0;
2164
2165	for(int ict=0;ict<ct_Number;ict++){
2166
2167	// Setting trigger conditions
2168	Trigger_CT[ict]->SetMultiplicity(Trigger_multiplicity[ict]);
2169	Trigger_CT[ict]->SetTopology(Trigger_topology[ict]);
2170	for (int i=0;i<ct_NPixels;i++)
2171	fpixelthres[i]=qThreshold[ict][i];
2172
2173	// Rise the discrimnator threshold to avoid huge rates
2174	if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
2175	for(int ii=0;ii<ct_NPixels;ii++){
2176	if( nsb_phepns_rotated[ict][ii]>riseDiskThres)
2177	fpixelthres[ii]=secureDiskThres;
2178	}
2179
2180	Trigger_CT[ict]->SetThreshold(fpixelthres);
2181
2182	Trigger_CT[ict]->Diskriminate() ;
2183
2184	//
2185	// look if in all the signals in the trigger signal branch
2186	// is a possible Trigger. Therefore we habe to diskriminate all
2187	// the simulated analog signals (Method Diskriminate in class
2188	// MTrigger). We look simultanously for the moments at which
2189	// there are more than TRIGGER_MULTI pixels above the
2190	// CHANNEL_THRESHOLD.
2191	//
2192
2193	Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ;
2194
2195	Lev1MT[ict] = 0 ;
2196
2197	//
2198	// Start the First Level Trigger simulation
2199	//
2200
2201	if ( Lev0MT[ict] > 0 \|\| Write_All_Images) {
2202
2203	Lev1MT[ict]= Trigger_CT[ict]->FirstLevel();
2204	}
2205	if (Lev1MT[ict]>0){
2206	++ntrigger[ict];
2207	}
2208	}
2209
2210	Int_t NumImages = 0;
2211	Int_t CT_triggered=0;
2212	for(int ict=0;ict<ct_Number;ict++){
2213	if(NumImages==0 && Lev1MT[ict]>0)
2214	CT_triggered=ict;
2215	NumImages = (NumImages>=Lev1MT[ict])?NumImages:1;
2216	Lev0MT[ict]=1;
2217	if (Lev1MT[ict]==0 && Write_All_Images){
2218	NumImages=1;
2219	Lev0MT[ict]=0;
2220	}
2221	}
2222
2223	for(int ict=0;ict<ct_Number;ict++){
2224	for(Int_t ii=0;ii<NumImages;ii++){
2225
2226	btrigger=1;
2227
2228	// Loop over different level one triggers
2229
2230	//
2231	// fill inside class fadc the member output
2232	//
2233	if(Lev1MT[ict]>0)
2234	Fadc_CT[ict]->
2235	TriggeredFadc(Trigger_CT[ict]->GetFirstLevelTime(ii));
2236	else
2237	Fadc_CT[ict]->
2238	TriggeredFadc(Trigger_CT[CT_triggered]->GetFirstLevelTime(ii));
2239
2240	if(!flagstoring)
2241	nstoredevents++;
2242	flagstoring = 1;
2243
2244	if (Write_McTrig){
2245	McTrig[ict]->SetFirstLevel (Lev1MT[ict]);
2246	McTrig[ict]
2247	->SetTime(Trigger_CT[ict]->GetFirstLevelTime(ii),ii+1);
2248	Trigger_CT[ict]->GetMapDiskriminator(trigger_map);
2249	McTrig[ict]->SetMapPixels(trigger_map,ii);
2250	}
2251
2252	// Fill Evt information
2253
2254	if (Write_RawEvt){
2255
2256	//
2257	// Fill the header of this event
2258	//
2259
2260	EvtHeader[ict]
2261	->FillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ;
2262
2263	// fill pixel information
2264
2265	if (Lev1MT[ict] \|\| Write_All_Images){
2266	if (addElecNoise) Fadc_CT[ict]->DigitalNoise();
2267	for(UInt_t i=0;
2268	i<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
2269	i++){
2270	if(!Fadc_CT[ict]->IsPixelUsed(i)) continue;
2271	for (j=0;j<FADC_SLICES;j++){
2272	fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(i,j),j);
2273	fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(i,j),j);
2274	}
2275	EvtData[ict]->AddPixel(i,fadcValues,0);
2276	EvtData[ict]->AddPixel(i,fadcValuesLow,kTRUE);
2277	}
2278	}
2279	}
2280	}
2281	//
2282	// if a first level trigger occurred, then
2283	// 1. do some other stuff (not implemented)
2284	// 2. start the gui tool
2285
2286	if(FADC_Scan){
2287	if ( Lev0MT[ict] > 0 ) {
2288	Fadc_CT[ict]->ShowSignal( McEvt[ict], (Float_t) 60. ) ;
2289	}
2290	}
2291
2292	if(Trigger_Scan){
2293	if ( Lev0MT[ict] > 0 ) {
2294	Trigger_CT[ict]->ShowSignal(McEvt[ict]) ;
2295	}
2296	}
2297
2298	} // end CT loop
2299
2300	// If there is trigger in some telecope or we store all showers
2301	if(btrigger){
2302	if (Write_McEvt){
2303	//
2304	// fill the MMcEvt with all information
2305	//
2306
2307	for (int ict=0;ict<ct_Number;ict++){
2308	Float_t ftime, ltime;
2309	if (reflector_file_version<6){
2310	mcevth[ict].get_times(&ftime, &ltime);
2311	McEvt[ict]->Fill( 0,
2312	(UShort_t) mcevth[0].get_primary() ,
2313	mcevth[ict].get_energy(),
2314	-1.0,
2315	-1.0,
2316	-1.0,
2317	mcevth[ict].get_theta(),
2318	mcevth[ict].get_phi(),
2319	mcevth[ict].get_core(),
2320	mcevth[ict].get_coreX(),
2321	mcevth[ict].get_coreY(),
2322	impactD[ict],
2323	phiCT[ict],
2324	thetaCT[ict],
2325	ftime,
2326	ltime,
2327	0,
2328	0,
2329	0,
2330	0,
2331	0,
2332	0,
2333	0,
2334	(UInt_t)mcevth[ict].get_CORSIKA(),
2335	(UInt_t)mcevth[ict].get_AtmAbs(),
2336	(UInt_t)(mcevth[ict].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()),
2337	(UInt_t) ncph[ict],
2338	(UInt_t) inumphe_CT[ict],
2339	(UInt_t) inumphensb[ict]+inumphe_CT[ict],
2340	-1.0,
2341	-1.0,
2342	-1.0);
2343	}
2344	else{
2345	Float_t Nmax, t0, tmax, a, b, c, chi2;
2346	mcevth_2[ict].get_times(&ftime, &ltime);
2347	chi2=mcevth_2[ict].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2348	McEvt[ict]->Fill( (UInt_t) mcevth_2[ict].get_evt_number(),
2349	(UShort_t) mcevth_2[ict].get_primary() ,
2350	mcevth_2[ict].get_energy(),
2351	mcevth_2[ict].get_thick0(),
2352	mcevth_2[ict].get_first_target(),
2353	mcevth_2[ict].get_z_first_int(),
2354	mcevth_2[ict].get_theta(),
2355	mcevth_2[ict].get_phi(),
2356	mcevth_2[ict].get_core(),
2357	mcevth_2[ict].get_coreX(),
2358	mcevth_2[ict].get_coreY(),
2359	impactD[ict],
2360	mcevth_2[ict].get_phi_CT(),
2361	mcevth_2[ict].get_theta_CT(),
2362	ftime,
2363	ltime,
2364	Nmax,
2365	t0,
2366	tmax,
2367	a,
2368	b,
2369	c,
2370	chi2,
2371	(UInt_t)mcevth_2[ict].get_CORSIKA(),
2372	(UInt_t)mcevth_2[ict].get_AtmAbs(),
2373	(UInt_t) (mcevth_2[ict].get_MirrAbs()+mcevth_2[ict].get_OutOfMirr()+mcevth_2[ict].get_BlackSpot()),
2374	(UInt_t) ncph[ict],
2375	(UInt_t) inumphe_CT[ict],
2376	(UInt_t) inumphensb[ict]+inumphe_CT[ict],
2377	mcevth_2[ict].get_ElecFraction(),
2378	mcevth_2[ict].get_MuonFraction(),
2379	mcevth_2[ict].get_OtherFraction());
2380	}
2381	}
2382	}
2383	// We don not count photons out of the camera.
2384
2385
2386	//
2387	// write it out to the file outfile
2388	//
2389
2390	EvtTree.Fill() ;
2391	}
2392
2393	// clear all
2394	for(int ict=0;ict<ct_Number;ict++){
2395	if (Write_RawEvt) EvtHeader[ict]->Clear() ;
2396	if (Write_RawEvt) EvtData[ict]->DeletePixels();
2397	if (Write_McTrig) McTrig[ict]->Clear() ;
2398	if (Write_McEvt) McEvt[ict]->Clear() ;
2399	}
2400	}
2401
2402	#ifdef __DEBUG__
2403	printf("\n");
2404
2405	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
2406
2407	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
2408
2409	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
2410
2411	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
2412
2413	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
2414	(int)pixels[ici][icj][PIXNUM],
2415	pixels[ici][icj][PIXX],
2416	pixels[ici][icj][PIXY],
2417	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
2418
2419	}
2420	}
2421	}
2422	}
2423
2424	for (int i=0;
2425	i<((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels(); ++i) {
2426	printf("%d (%d): ", i, npixneig[i]);
2427	for (j=0; j<npixneig[i]; ++i)
2428	printf(" %d", pixneig[i][j]);
2429	printf("\n");
2430	}
2431
2432	#endif // __DEBUG__
2433
2434
2435	// We search the maximum impact parameter fo the simualted showers
2436	maxpimpact=maxpimpact<impactD[0]?impactD[0]:maxpimpact;
2437
2438	// look for the next event
2439
2440	for(int ict=0;ict<ct_Number;ict++)
2441	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2442
2443	} // end while there is a next event
2444
2445	if( !isA( flag, FLAG_END_OF_RUN )){
2446	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
2447	break;
2448	}
2449	else { // found end of run
2450	ntshow += nshow;
2451	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
2452
2453	//fread( flag, SIZE_OF_FLAGS, 1, inputfile );
2454	for(int ict=0;ict<ct_Number;ict++)
2455	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2456
2457	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
2458	log(SIGNATURE, "End of file . . .\n");
2459	still_in_loop = FALSE;
2460	log(SIGNATURE, "Reading ASCII files at the end of the reflector file. . .\n");
2461	for(int ict=0;ict<ct_Number;ict++){
2462	read_ascii(inputfile[ict], McConfigRunHeader[ict]);
2463	}
2464	if ((! Data_From_STDIN) && ( !feof(inputfile[0]) )){
2465
2466	// we have concatenated input files.
2467	// get signature of the next part and check it.
2468
2469	if((reflector_file_version=check_reflector_file( inputfile[0] ))==FALSE){
2470	log(SIGNATURE, "Next file is not recognised as a reflector file.\n");
2471	log(SIGNATURE, "Stopping ...\n");
2472	break;
2473	}
2474
2475	}
2476
2477	for(int ict=0;ict<ct_Number;ict++)
2478	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2479
2480	} // end if found end of file
2481	} // end if found end of run
2482
2483	} // end if else found start of run
2484	} // end big while loop
2485
2486	//<@ Finally we should fill th McRunHeader
2487
2488	Float_t heights[10];
2489	time_t ltime;
2490	Float_t ftime;
2491	Float_t rnum;
2492	Float_t viewcone[2]={0,0};
2493
2494	get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
2495	if (reflector_file_version<6){
2496	telestheta=-10.0;
2497	telesphi=-10.0;
2498	mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2499	mcevth[0].get_phi_range(&shphimin,&shphimax);
2500	mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2501	mcevth[0].get_phi_range(&shphimin,&shphimax);
2502	corsika=UInt_t(mcevth[0].get_VersionPGM()*1000);
2503	for (int i=0; i< 10;i++)
2504	heights[i]=mcevth[0].get_HeightLev (i);
2505	rnum=mcevth[0].get_RunNumber();
2506	}
2507	else{
2508	telestheta=mcevth_2[0].get_theta_CT();
2509	telesphi=mcevth_2[0].get_phi_CT();
2510	mcevth_2[0].get_theta_range(&shthetamin, &shthetamax);
2511	mcevth_2[0].get_phi_range(&shphimin,&shphimax);
2512	corsika=UInt_t(mcevth_2[0].get_VersionPGM()*1000);
2513	for (int i=0; i< 10;i++)
2514	heights[i]=mcevth_2[0].get_HeightLev (i);
2515	rnum=mcevth_2[0].get_RunNumber();
2516	mcevth_2[0].get_viewcone(&viewcone[0],&viewcone[1]);
2517	}
2518
2519	get_source_off(&sofftheta,&soffphi);
2520	if(!Trigger_Loop) icontrigger=0;
2521	time (&ltime);
2522	ftime = ((Float_t)ltime)/1000;
2523
2524	if (reflector_file_version<6)
2525	McRunHeader->Fill(rnum,
2526	(UInt_t) 0,
2527	mcevth[0].get_DateRun(),
2528	ftime,
2529	icontrigger,
2530	!Write_All_Images,
2531	Write_McEvt,
2532	Write_McTrig,
2533	Write_McFADC,
2534	Write_RawEvt,
2535	addElecNoise,
2536	ct_NPixels,
2537	(UInt_t)ntshow,
2538	(UInt_t)nstoredevents,
2539	0,
2540	sfRaH,
2541	sfRaM,
2542	sfRaS,
2543	sfDeD,
2544	sfDeM,
2545	sfDeS,
2546	meanNSB,
2547	telestheta,
2548	telesphi,
2549	sofftheta,
2550	soffphi,
2551	shthetamax,
2552	shthetamin,
2553	shphimax,
2554	shphimin,
2555	maxpimpact,
2556	mcevth[0].get_CWaveLower(),
2557	mcevth[0].get_CWaveUpper(),
2558	mcevth[0].get_slope(),
2559	1,
2560	heights,
2561	corsika,
2562	(UInt_t)(reflector_file_version*100),
2563	(UInt_t)(VERSION*100),
2564	0);
2565	else
2566	McRunHeader->Fill(rnum,
2567	(UInt_t) 0,
2568	mcevth_2[0].get_DateRun(),
2569	ftime,
2570	icontrigger,
2571	!Write_All_Images,
2572	Write_McEvt,
2573	Write_McTrig,
2574	Write_McFADC,
2575	Write_RawEvt,
2576	addElecNoise,
2577	ct_NPixels,
2578	(UInt_t)ntshow,
2579	(UInt_t)nstoredevents,
2580	0,
2581	sfRaH,
2582	sfRaM,
2583	sfRaS,
2584	sfDeD,
2585	sfDeM,
2586	sfDeS,
2587	meanNSB,
2588	telestheta,
2589	telesphi,
2590	sofftheta,
2591	soffphi,
2592	shthetamax,
2593	shthetamin,
2594	shphimax,
2595	shphimin,
2596	maxpimpact,
2597	mcevth_2[0].get_CWaveLower(),
2598	mcevth_2[0].get_CWaveUpper(),
2599	mcevth_2[0].get_slope(),
2600	1,
2601	heights,
2602	corsika,
2603	(UInt_t)(reflector_file_version*100),
2604	(UInt_t)(VERSION*100),
2605	0);
2606	// Fill some missing values for MRawRunHeader
2607
2608	RunHeader->SetRunNumber(rnum);
2609	RunHeader->SetNumEvents(nstoredevents);
2610
2611	// Fill MMcCorsikaRunHeader
2612
2613	Float_t constantC[50];
2614	mcrunh.get_constantC(&constantC[0]);
2615	Float_t constantCKA[40];
2616	mcrunh.get_constantCKA(&constantCKA[0]);
2617	Float_t constantCETA[5];
2618	mcrunh.get_constantCETA(&constantCETA[0]);
2619	Float_t constantCSTRBA[11];
2620	mcrunh.get_constantCSTRBA(&constantCSTRBA[0]);
2621	Float_t constantAATM[5];
2622	mcrunh.get_constantAATM(&constantAATM[0]);
2623	Float_t constantBATM[5];
2624	mcrunh.get_constantBATM(&constantBATM[0]);
2625	Float_t constantCATM[5];
2626	mcrunh.get_constantCATM(&constantCATM[0]);
2627	Float_t constantNFL[4];
2628	mcrunh.get_constantNFL(&constantNFL[0]);
2629
2630	if(reflector_file_version>5)
2631	McCorsikaRunHeader->Fill(rnum,
2632	mcrunh.get_date(),
2633	corsika,
2634	1,
2635	heights,
2636	mcevth_2[0].get_slope(),
2637	mcrunh.get_ELow(),
2638	mcrunh.get_EUpp(),
2639	mcrunh.get_EGS4(),
2640	mcrunh.get_NKG(),
2641	mcrunh.get_Ecutoffh(),
2642	mcrunh.get_Ecutoffm(),
2643	mcrunh.get_Ecutoffe(),
2644	mcrunh.get_Ecutoffg(),
2645	constantC,
2646	constantCKA,
2647	constantCETA,
2648	constantCSTRBA,
2649	constantAATM,
2650	constantBATM,
2651	constantCATM,
2652	constantNFL,
2653	viewcone,
2654	mcrunh.get_wobble(),
2655	mcrunh.get_atmophere()
2656	);
2657
2658	// Store qe for each PMT in output file
2659	TArrayF qe_pmt;
2660	TArrayF wav_pmt;
2661
2662	for(int ict=0;ict<ct_Number;ict++){
2663	McConfigRunHeader[ict]->InitSizePMTs(ct_NPixels);
2664	for(int i=0; i<(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();i++){
2665	McConfigRunHeader[ict]->AddPMT(i);
2666	MGeomPMT &pmt = McConfigRunHeader[ict]->GetPMT(i);
2667	qe_pmt.Set(pointsQE[ict],QE[ict][i][1]);
2668	wav_pmt.Set(pointsQE[ict],QElambda);
2669	pmt.SetArraySize(pointsQE[ict]);
2670	pmt.SetPMTContent(i,wav_pmt,qe_pmt);
2671	}
2672
2673	// Store Light Guides factors in the output file
2674	TArrayF theta_lg;
2675	TArrayF factor_lg;
2676	theta_lg.Set(pointsWC,WC[0]);
2677	factor_lg.Set(pointsWC,WC[1]);
2678	McConfigRunHeader[ict]->SetLightGuides(theta_lg,factor_lg);
2679
2680	}
2681
2682	// Fill the Header Tree with the current leaves of each branch
2683	HeaderTree.Fill() ;
2684
2685	//++
2686	// put the Event to the root file
2687	//--
2688
2689	outfile_temp.Write() ;
2690	outfile_temp.Close() ;
2691
2692	// close input file
2693
2694	if (Trigger_Loop){
2695	log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
2696	ntshow, ntcph[0] );
2697	datafile<<ntshow<<" event(s), with a total of "<<ntcph[0]<<" C.photons"<<endl;
2698	log( SIGNATURE, "Trigger Mode. \n");
2699	log( SIGNATURE, "Fraction of triggers: \n");
2700	datafile<<"Fraction of triggers: "<<endl;
2701	for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
2702	for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
2703	for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
2704	log( SIGNATURE, "Thres %5.1f, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
2705	fthrescount,imulticount,isorttopo[itopocount],((float)ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0), ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop], ntshow);
2706	datafile<<"Thres "<<fthrescount<<", Multi "<<imulticount<<", Topo"<<isorttopo[itopocount]<<": ";
2707	datafile<<((float)ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0)<<"% ("<<ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]<<" out of "<<ntshow<<")"<<endl;
2708	}
2709	}
2710	}
2711	}
2712	else{
2713	for(int ict=0;ict<ct_Number;ict++){
2714	log( SIGNATURE,
2715	"%d event(s), with a total of %d C.photons in CT %i (%s)\n",
2716	ntshow, ntcph[ict],ict,GeometryName[ict] );
2717	datafile<<ntshow<<" event(s), with a total of "<<ntcph[ict]
2718	<<" C.photons in CT "<<ict<<" ("<<GeometryName[ict]<<")"<<endl;
2719	log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
2720	((float)ntrigger[ict]) / ((float)ntshow) * 100.0, ntrigger[ict], ntshow);
2721	datafile<<"Fraction of triggers: "<<((float)ntrigger[ict]) / ((float)ntshow) * 100.0<<" ("<<ntrigger[ict]<<" out of "<<ntshow<<" )"<<endl;
2722	}
2723	}
2724
2725	// close files
2726
2727	log( SIGNATURE, "Closing files\n" );
2728
2729	if( ! Data_From_STDIN ){
2730	for(int ict=0;ict<ct_Number;ict++)
2731	fclose( inputfile[ict] );
2732	}
2733	datafile.close();
2734
2735	// program finished
2736
2737	log( SIGNATURE, "Done.\n");
2738
2739	return( 0 );
2740	}
2741	//!@}
2742
2743	// @T \newpage
2744
2745	//!@subsection Functions definition.
2746
2747	//!-----------------------------------------------------------
2748	// @name present
2749	//
2750	// @desc Make some presentation
2751	//
2752	// @date Sat Jun 27 05:58:56 MET DST 1998
2753	//------------------------------------------------------------
2754	// @function
2755
2756	//!@{
2757	void
2758	present(void)
2759	{
2760	cout << "##################################################\n"
2761	<< SIGNATURE << '\n' << '\n'
2762	<< "Processor of the reflector output\n"
2763	<< "J C Gonzalez, Jun 1998\n"
2764	<< "##################################################\n\n"
2765	<< flush ;
2766	}
2767	//!@}
2768
2769
2770	//!-----------------------------------------------------------
2771	// @name usage
2772	//
2773	// @desc show help
2774	//
2775	// @date Tue Dec 15 16:23:30 MET 1998
2776	//------------------------------------------------------------
2777	// @function
2778
2779	//!@{
2780	void
2781	usage(void)
2782	{
2783	present();
2784	cout << "\nusage ::\n\n"
2785	<< "\t camera "
2786	<< " [ -@ paramfile ] "
2787	<< " [ -h ] "
2788	<< "\n\n or \n\n"
2789	<< "\t camera < paramfile"
2790	<< "\n\n";
2791	exit(0);
2792	}
2793	//!@}
2794
2795
2796	//!-----------------------------------------------------------
2797	// @name log
2798	//
2799	// @desc function to send log information
2800	//
2801	// @var funct Name of the caller function
2802	// @var fmt Format to be used (message)
2803	// @var ... Other information to be shown
2804	//
2805	// @date Sat Jun 27 05:58:56 MET DST 1998
2806	//------------------------------------------------------------
2807	// @function
2808
2809	//!@{
2810	void
2811	log(const char funct, char fmt, ...)
2812	{
2813	va_list args;
2814
2815	// Display the name of the function that called error
2816	printf("[%s]: ", funct);
2817
2818	// Display the remainder of the message
2819	va_start(args, fmt);
2820	vprintf(fmt, args);
2821	va_end(args);
2822	}
2823	//!@}
2824
2825
2826	//!-----------------------------------------------------------
2827	// @name error
2828	//
2829	// @desc function to send an error message, and abort the program
2830	//
2831	// @var funct Name of the caller function
2832	// @var fmt Format to be used (message)
2833	// @var ... Other information to be shown
2834	//
2835	// @date Sat Jun 27 05:58:56 MET DST 1998
2836	//------------------------------------------------------------
2837	// @function
2838
2839	//!@{
2840	void
2841	error(const char funct, char fmt, ...)
2842	{
2843	va_list args;
2844
2845	// Display the name of the function that called error
2846	fprintf(stdout, "ERROR in %s: ", funct);
2847
2848	// Display the remainder of the message
2849	va_start(args, fmt);
2850	vfprintf(stdout, fmt, args);
2851	va_end(args);
2852
2853	perror(funct);
2854
2855	exit(1);
2856	}
2857	//!@}
2858
2859
2860	//!-----------------------------------------------------------
2861	// @name isA
2862	//
2863	// @desc returns TRUE(FALSE), if the flag is(is not) the given
2864	//
2865	// @var s1 String to be searched
2866	// @var flag Flag to compare with string s1
2867	// @return TRUE: both strings match; FALSE: oth.
2868	//
2869	// @date Wed Jul 8 15:25:39 MET DST 1998
2870	//------------------------------------------------------------
2871	// @function
2872
2873	//!@{
2874	int
2875	isA( char * s1, const char * flag ) {
2876	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
2877	}
2878	//!@}
2879
2880
2881	//!-----------------------------------------------------------
2882	// @name read_QE
2883	//
2884	// @desc read QE data
2885	//
2886	// @date thu 5 17:59:57 CEST 2002
2887	//------------------------------------------------------------
2888	// @function
2889
2890	//!@{
2891	void
2892	read_QE(char fname[256], int ict){
2893	ifstream qefile;
2894	char line[LINE_MAX_LENGTH];
2895	int i, j, icount;
2896	float qe;
2897
2898	//------------------------------------------------------------
2899	// second, pixels' QE
2900
2901	// try to open the file
2902
2903	log("read_QE", "Opening the file \"%s\" . . .\n", fname);
2904
2905	qefile.open( fname );
2906
2907	// if it is wrong or does not exist, exit
2908
2909	if ( qefile.bad() )
2910	error( "read_QE", "Cannot open \"%s\". Exiting.\n", fname );
2911
2912	// read file
2913
2914	log("read_QE", "Reading data . . .\n");
2915
2916	i=-1;
2917	icount = 0;
2918
2919	while ( ! qefile.eof() ) {
2920
2921	// get line from the file
2922
2923	qefile.getline(line, LINE_MAX_LENGTH);
2924
2925	// skip if comment
2926
2927	if ( *line == '#' )
2928	continue;
2929
2930	// if it is the first valid value, it is the number of QE data points
2931
2932	if ( i < 0 ) {
2933
2934	// get the number of datapoints
2935
2936	sscanf(line, "%d", &pointsQE[ict]);
2937
2938	// allocate memory for the table of QEs
2939
2940	QE[ict] = new float ** [ct_NPixels];
2941
2942	for ( i=0; i<ct_NPixels; ++i ) {
2943	QE[ict][i] = new float * [2];
2944	QE[ict][i][0] = new float[pointsQE[ict]];
2945	QE[ict][i][1] = new float[pointsQE[ict]];
2946	}
2947
2948	QElambda = new float [pointsQE[ict]];
2949
2950	for ( i=0; i<pointsQE[ict]; ++i ) {
2951	qefile.getline(line, LINE_MAX_LENGTH);
2952	sscanf(line, "%f", &QElambda[i]);
2953	}
2954
2955	i=0;
2956
2957	continue;
2958	}
2959
2960	// get the values (num-pixel, num-datapoint, QE-value)
2961
2962	if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2963	break;
2964
2965	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2966	((j-1) < pointsQE[ict]) && ((j-1) > -1) ) {
2967	QE[ict][i-1][0][j-1] = QElambda[j-1];
2968	QE[ict][i-1][1][j-1] = qe;
2969	}
2970
2971	if ( i > ct_NPixels) break;
2972
2973	icount++;
2974
2975	}
2976
2977	if(icount/pointsQE[ict] < ct_NPixels){
2978	error( "read_QE", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2979	icount/pointsQE[ict], ct_NPixels );
2980	}
2981
2982	// close file
2983
2984	qefile.close();
2985
2986	// test QE
2987
2988	for(icount=0; icount< ct_NPixels; icount++){
2989	for(i=0; i<pointsQE[ict]; i++){
2990	if( QE[ict][icount][0][i] < 100. \|\| QE[ict][icount][0][i] > 1000. \|\|
2991	QE[ict][icount][1][i] < 0. \|\| QE[ict][icount][1][i] > 100.){
2992	error( "read_QE", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2993	icount, i, QE[ict][icount][0][i], QE[ict][icount][1][i] );
2994	}
2995	}
2996	}
2997
2998	// end
2999
3000	log("read_QE", "Done.\n");
3001	}
3002	//!@}
3003
3004	//!-----------------------------------------------------------
3005	// @name read_WC
3006	//
3007	// @desc read WC data
3008	//
3009	// @date thu 5 17:59:57 CEST 2002
3010	//------------------------------------------------------------
3011	// @function
3012
3013	//!@{
3014	void
3015	read_WC(void){
3016	ifstream wcfile;
3017	char line[LINE_MAX_LENGTH];
3018	int i;
3019
3020	//------------------------------------------------------------
3021	// Read Light Guides data
3022
3023	// try to open the file
3024
3025	log("read_WC", "Opening the file \"%s\" . . .\n", WC_FILE);
3026
3027	wcfile.open( WC_FILE );
3028
3029	// if it is wrong or does not exist, exit
3030
3031	if ( wcfile.bad() )
3032	error( "read_WC", "Cannot open \"%s\". Exiting.\n", WC_FILE );
3033
3034	// read file
3035
3036	log("read_WC", "Reading data . . .\n");
3037
3038	// get line from the file
3039
3040	wcfile.getline(line, LINE_MAX_LENGTH);
3041
3042	// get the number of datapoints
3043
3044	sscanf(line, "%d", &pointsWC);
3045
3046	// allocate memory for the table of QEs
3047
3048	WC = new float * [2];
3049	WC[0] = new float[pointsWC];
3050	WC[1] = new float[pointsWC];
3051
3052	for ( i=0; i<pointsWC; ++i ) {
3053	wcfile.getline(line, LINE_MAX_LENGTH);
3054	sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
3055	}
3056
3057	// close file
3058
3059	wcfile.close();
3060
3061	// read
3062
3063	log("read_WC", "Done.\n");
3064	}
3065	//!@}
3066
3067
3068	//!-----------------------------------------------------------
3069	// @name read_ascii
3070	//
3071	// @desc read ascii configuration files used by the reflector
3072	//
3073	// @date tue dec 10 17:14:10 CET 2002
3074	//------------------------------------------------------------
3075	// @function
3076
3077	//!@{
3078	void
3079	read_ascii(FILE sp, MMcConfigRunHeader config)
3080	{
3081	Float_t radius = -1.0;
3082	Float_t focal = -1.0;
3083	Float_t stdfocal = -1.0;
3084	Float_t point = -1.0;
3085	Float_t stdpoint = -1.0;
3086	Float_t adjust = -1.0;
3087	Float_t spot = -1.0;
3088	Float_t camwidth = -1.0;
3089
3090	Int_t imir;
3091	Float_t f,sx,sy,x,y,z,thetan,phin,xn,yn,zn;
3092	Float_t dx,dy;
3093	Int_t nummir, numref;
3094	Float_t wav,ref;
3095
3096	Char_t token[40];
3097	Char_t line[511];
3098	Char_t flag;
3099
3100	while(1){
3101	if((flag=fgetc(sp))==EOF)
3102	break;
3103	fgets(&line[1],500,sp);
3104	line[0]=flag;
3105	if (line[0]== '#' && line[2]=='n' && line[3]=='u' && line[4]=='m' &&
3106	line[5]== 'b' && line[12]=='d' && line[13]=='a' && line[14]=='t'){
3107	fgets(line,500,sp);
3108	sscanf(line,"%i",&numref);
3109	TArrayF wavarray(numref);
3110	TArrayF refarray(numref);
3111	fgets(line,500,sp);
3112	for(int i=0; i<numref;i++){
3113	fgets(line,500,sp);
3114	sscanf(line,"%f %f",&wav,&ref);
3115	wavarray[i]=wav;
3116	refarray[i]=ref;
3117	}
3118	for (int j=0; j<nummir;j++){
3119	MGeomMirror &mirror = config->GetMirror(j);
3120	mirror.SetArraySize(numref);
3121	mirror.SetReflectivity(wavarray, refarray);
3122	}
3123	continue;
3124	}
3125	if (line[0]== '#')
3126	continue;
3127	if (strstr(line, "type") == line)
3128	continue;
3129	if (strstr(line, "focal_distance") == line){
3130	sscanf(line,"%s %f",token,&focal);
3131	continue;
3132	}
3133	if (strstr(line, "focal_std") == line){
3134	sscanf(line,"%s %f",token,&stdfocal);
3135	continue;
3136	}
3137	if (strstr(line, "point_spread") == line){
3138	sscanf(line,"%s %f",token,&point);
3139	continue;
3140	}
3141	if (strstr(line, "point_std") == line){
3142	sscanf(line,"%s %f",token,&stdpoint);
3143	continue;
3144	}
3145	if (strstr(line, "adjustment_dev") == line){
3146	sscanf(line,"%s %f",token,&adjust);
3147	continue;
3148	}
3149	if (strstr(line, "black_spot") == line){
3150	sscanf(line,"%s %f",token,&spot);
3151	continue;
3152	}
3153	if (strstr(line, "camera_width") == line){
3154	sscanf(line,"%s %f",token,&camwidth);
3155	continue;
3156	}
3157	if (strstr(line, "n_pixels") == line)
3158	continue;
3159	if (strstr(line, "pixel_width") == line)
3160	continue;
3161	if (strstr(line, "n_centralpixels") == line)
3162	continue;
3163	if (strstr(line, "n_gappixels") == line)
3164	continue;
3165	if (strstr(line, "n_mirrors") == line){
3166	sscanf(line,"%s %i",token,&nummir);
3167	config->InitSizeMirror(nummir);
3168	continue;
3169	}
3170	if (strstr(line, "r_mirror") == line){
3171	sscanf(line,"%s %f",token,&radius);
3172	continue;
3173	}
3174	if (strstr(line, "define_mirrors") == line){
3175	for(int i=0;i<nummir;i++){
3176	fgets(line,500,sp);
3177	sscanf(line,"%i %f %f %f %f %f %f %f %f %f %f %f",
3178	&imir,&f,&sx,&sy,&x,&y,&z,&thetan,&phin,&xn,&yn,&zn);
3179	config->AddMirror(i);
3180	MGeomMirror &mirror = config->GetMirror(i);
3181	mirror.SetMirrorContent(imir,f,sx,sy,x,y,z,thetan,phin,xn,yn,zn);
3182	}
3183	fgets(line,500,sp);
3184	while(line[2]!='w' \|\| line[3]!='i' \|\| line[4]!='t' \|\| line[5]!='h'){
3185	fgets(line,500,sp);
3186	}
3187	fgets(line,500,sp);
3188	for(int i=0;i<nummir;i++){
3189	fgets(line,500,sp);
3190	sscanf(line,"%f %f",&dx,&dy);
3191	MGeomMirror &mirror = config->GetMirror(i);
3192	mirror.SetMirrorDeviations(dx,dy);
3193	}
3194	continue;
3195	}
3196	}
3197	config->SetMagicDef(radius, focal, stdfocal, point,
3198	stdpoint, adjust, spot, camwidth);
3199	}
3200
3201
3202	//!-----------------------------------------------------------
3203	// @name igen_pixel_coordinates
3204	//
3205	// @desc generate the pixel center coordinates
3206	//
3207	// @var *pcam structure camera containing all the
3208	// camera information
3209	// @return total number of pixels
3210	//
3211	// DP
3212	//
3213	// @date Thu Oct 14 10:41:03 CEST 1999
3214	//------------------------------------------------------------
3215	// @function
3216
3217	//!@{
3218	/****** igen_pixel_coordinates() *******************************/
3219
3220	int igen_pixel_coordinates(struct camera *pcam) {
3221	/* generate pixel coordinates, return value is number of pixels */
3222
3223	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
3224	float fsegment_fract;
3225	double dtsize;
3226	double dhsize;
3227	double dpsize;
3228	double dxfirst_pix;
3229	double dyfirst_pix;
3230	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5;
3231	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5;
3232
3233
3234	double dstartx, dstarty; /* for the gap pixels and outer pixels */
3235	int j, nrow;
3236
3237	dpsize = pcam->dpixdiameter_cm;
3238	dtsize = dpsize * sqrt(3.) / 2.;
3239	dhsize = dpsize / 2.;
3240
3241	/* Loop over central pixels to generate co-ordinates */
3242
3243	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
3244
3245	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
3246
3247	pcam->dpixsizefactor[ipixno-1] = 1.;
3248
3249	in = 0;
3250
3251	i = 0;
3252	itot_inside_ring = 0;
3253	iring_no = 0;
3254
3255	/* Calculate the number of pixels out to and including the ring containing pixel number */
3256	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
3257
3258	while (itot_inside_ring == 0){
3259
3260	iN = 3(i(i+1)) + 1;
3261
3262	if (ipixno <= iN){
3263	iring_no = i;
3264	itot_inside_ring = iN;
3265	}
3266
3267	i++;
3268	}
3269
3270
3271	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
3272
3273	ipix_in_ring = 0;
3274	for (i = 0; i < iring_no; ++i){
3275
3276	ipix_in_ring = ipix_in_ring + 6;
3277	}
3278
3279	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
3280	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
3281	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
3282
3283	isegment = 0;
3284	fsegment_fract = 0.;
3285	if (iring_no > 0) {
3286
3287	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
3288
3289	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
3290
3291	}
3292
3293	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
3294	/* pixel width (flat to flat) dpsize. */
3295
3296	dxfirst_pix = dpsize*iring_no;
3297	dyfirst_pix = 0.;
3298
3299	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
3300
3301	ddxseg1 = - dhsize*iring_no;
3302	ddyseg1 = dtsize*iring_no;
3303	ddxseg2 = -dpsize*iring_no;
3304	ddyseg2 = 0.;
3305	ddxseg3 = ddxseg1;
3306	ddyseg3 = -ddyseg1;
3307	ddxseg4 = -ddxseg1;
3308	ddyseg4 = -ddyseg1;
3309	ddxseg5 = -ddxseg2;
3310	ddyseg5 = 0.;
3311
3312	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
3313	/* anti-clockwise around the ring by adding the segment to segment vectors. */
3314
3315	switch (isegment) {
3316
3317	case 0:
3318
3319	pcam->dxc[ipixno-1] = 0.;
3320	pcam->dyc[ipixno-1] = 0.;
3321
3322	case 1:
3323	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
3324	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
3325
3326	break;
3327
3328	case 2:
3329
3330	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
3331	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
3332
3333	break;
3334
3335	case 3:
3336
3337	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
3338	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
3339
3340	break;
3341
3342	case 4:
3343
3344	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
3345	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
3346
3347	break;
3348
3349	case 5:
3350
3351	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
3352	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
3353
3354	break;
3355
3356	case 6:
3357
3358	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
3359	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
3360
3361	break;
3362
3363	default:
3364
3365	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
3366	return(0);
3367
3368	} /* end switch */
3369
3370	} /* end for */
3371
3372	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
3373	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
3374
3375	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
3376
3377	j = pcam->inumcentralpixels;
3378
3379	for(i=0; i<pcam->inumgappixels; i=i+6){
3380	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
3381	pcam->dyc[j + i ] = dstarty;
3382	pcam->dpixsizefactor[j + i] = 1.;
3383	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
3384	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
3385	pcam->dpixsizefactor[j + i + 1] = 1.;
3386	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
3387	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
3388	pcam->dpixsizefactor[j + i+ 2] = 1.;
3389	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
3390	pcam->dyc[j + i + 3] = dstarty;
3391	pcam->dpixsizefactor[j + i+ 3] = 1.;
3392	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
3393	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
3394	pcam->dpixsizefactor[j + i+ 4] = 1.;
3395	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
3396	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
3397	pcam->dpixsizefactor[j + i + 5] = 1.;
3398	} /* end for */
3399	} /* end if */
3400
3401	/* generate positions of the outer pixels */
3402
3403	if( pcam->inumbigpixels > 0 ){
3404
3405	j = pcam->inumcentralpixels + pcam->inumgappixels;
3406
3407	for(i=0; i<pcam->inumbigpixels; i++){
3408	pcam->dpixsizefactor[j + i] = 2.;
3409	}
3410
3411	in = 0;
3412
3413	nrow = (int) ceil(dstartx / 2. / dpsize);
3414
3415	while(in < pcam->inumbigpixels){
3416
3417	pcam->dxc[j + in] = dstartx + dpsize;
3418	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
3419	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
3420	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
3421	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
3422	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
3423	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
3424	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
3425	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
3426	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
3427	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
3428	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
3429	for(i=1; i<nrow; i++){
3430	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
3431	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
3432	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
3433	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
3434	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
3435	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
3436	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
3437	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
3438	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
3439	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
3440	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
3441	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
3442	}
3443	in = in + 6 * nrow;
3444	dstartx = dstartx + 2. * dpsize;
3445	nrow = nrow + 1;
3446	} /* end while */
3447
3448	} /* end if */
3449
3450	return(pcam->inumpixels);
3451
3452	}
3453	//!@}
3454
3455	//!-----------------------------------------------------------
3456	// @name bpoint_is_in_pix
3457	//
3458	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
3459	//
3460	// @var *pcam structure camera containing all the
3461	// camera information
3462	// @var dx, dy point coordinates in centimeters
3463	// @var ipixnum pixel number (starting at 0)
3464	// @return TRUE if the point is inside the pixel, FALSE otherwise
3465	//
3466	// DP
3467	//
3468	// @date Thu Oct 14 16:59:04 CEST 1999
3469	//------------------------------------------------------------
3470	// @function
3471
3472	//!@{
3473
3474	/****** bpoint_is_in_pix() *************************************/
3475
3476	#define sqrt13 0.577350269 // = 1./sqrt(3.)
3477	#define sqrt3 1.732050807 // = sqrt(3.)
3478
3479	int bpoint_is_in_pix(double dx, double dy, MGeomCam *pgeomcam)
3480	{
3481	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
3482	/* the pixel is assumed to be a "closed set" */
3483
3484	/*
3485	a = length of one of the edges of one pixel,
3486	b = half the width of one pixel
3487	*/
3488
3489	const int numN = pgeomcam->GetNumPixels();
3490
3491	for (int i=0; i<numN; i++)
3492	{
3493	MGeomPix &pixel = (*pgeomcam)[i];
3494	const double b = pixel.GetD()/2;
3495	const double a = pixel.GetD()/sqrt3;
3496
3497	const double xx = dx - pixel.GetX();
3498	const double yy = dy - pixel.GetY();
3499
3500	if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) \|\|
3501	((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){
3502
3503	return i; // inside i
3504	}
3505
3506	// return -1; // outside
3507	}
3508
3509	return -1; // outside
3510	}
3511
3512	//!@}
3513
3514	//------------------------------------------------------------
3515	// @name dist_r_P
3516	//
3517	// @desc distance straight line r - point P
3518	//
3519	// @date Sat Jun 27 05:58:56 MET DST 1998
3520	// @function @code
3521	//------------------------------------------------------------
3522	// dist_r_P
3523	//
3524	// distance straight line r - point P
3525	//------------------------------------------------------------
3526
3527	float
3528	dist_r_P(float a, float b, float c,
3529	float l, float m, float n,
3530	float x, float y, float z)
3531	{
3532	return (
3533	sqrt((SQR((a-x)m-(b-y)l) +
3534	SQR((b-y)n-(c-z)m) +
3535	SQR((c-z)l-(a-x)n))/
3536	(SQR(l)+SQR(m)+SQR(n))
3537	)
3538	);
3539	}
3540
3541	//------------------------------------------------------------
3542	// @name check_reflector_file
3543	//
3544	// @desc check if a given reflector file has the right signature
3545	// @desc return TRUE or FALSE
3546	//
3547	// @date Mon Feb 14 16:44:21 CET 2000
3548	// @function @code
3549	//------------------------------------------------------------
3550
3551	int check_reflector_file(FILE *infile){
3552
3553	char sign[20]; // auxiliary variable
3554
3555	strcpy(sign, REFL_SIGNATURE_B);
3556
3557	fread( (char *)sign, strlen(REFL_SIGNATURE_B), 1, infile);
3558	if (strcmp(sign, REFL_SIGNATURE_A) == 0){
3559	fread( (char *)sign, 1, 1, infile);
3560	return 4;
3561	}
3562	else if (strcmp(sign, REFL_SIGNATURE_B) == 0){
3563	fread( (char *)sign, 1, 1, infile);
3564	return 5;
3565	}
3566	else if (strcmp(sign, REFL_SIGNATURE_C) == 0){
3567	// An empty bin has been removed and therefore we do not need to rezd it.
3568	return 6;
3569	}
3570	else {
3571	cout << "ERROR: Signature of .rfl file is not correct\n";
3572	cout << '"' << sign << '"' << '\n';
3573	cout << "should be: " << REFL_SIGNATURE_A <<" or "<< REFL_SIGNATURE_B <<" or " << REFL_SIGNATURE_C <<" or "<< '\n';
3574	return(FALSE);
3575	}
3576
3577
3578	}
3579
3580	//------------------------------------------------------------
3581	// @name lin_interpol
3582	//
3583	// @desc interpolate linearly between two points returning the
3584	// @desc y-value of the result
3585	//
3586	// @date Thu Feb 17 11:31:32 CET 2000
3587	// @function @code
3588	//------------------------------------------------------------
3589
3590	float lin_interpol( float x1, float y1, float x2, float y2, float x){
3591
3592	if( (x2 - x1)<=0. ){ // avoid division by zero, return average
3593	cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
3594	return((y1+y2)/2.);
3595	}
3596	else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
3597	return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
3598	}
3599	}
3600
3601	//------------------------------------------------------------
3602	// @name size_rotated
3603	//
3604	// @desc It rotates the NSB
3605	//
3606	// @date Tue Jan 7 17:09:25 CET 2003
3607	// @function @code
3608	//------------------------------------------------------------
3609
3610	int size_rotated(
3611	float *rotated,
3612	float nsb[iMAXNUMPIX],
3613	float rho)
3614	{
3615	int r=0;
3616	float size_rotated;
3617	float Num_Pixels;
3618	float PixNum=0;
3619	float rho_pixel;
3620	int j=0;
3621	int k=0;
3622
3623	for(int i=1; i<iMAXNUMPIX;i++){
3624	// Substituir per Int_t radius[iMAXNUMPIX]={0,1,1,1,1,1,1,2,2,2,2, ...}
3625	Num_Pixels=0;
3626	for (int ii=1; ii<17;ii++){
3627	if (Num_Pixels >= i){
3628	r=ii-1;
3629	break;
3630	}
3631
3632
3633
3634	if (ii<12){
3635	Num_Pixels+=ii*6;
3636	PixNum=6*ii;
3637	// size_rotated = (360/PixNum);
3638	//rho_pixel = rho/size_rotated;
3639	}
3640
3641
3642	else
3643	{
3644	Num_Pixels+=(ii-6)*6;
3645	PixNum=(ii-6)*6;
3646	// size_rotated = (360/PixNum);
3647	// rho_pixel = rho/size_rotated;
3648	}
3649	}
3650
3651
3652
3653	size_rotated = (360/PixNum);
3654	rho_pixel = rho/size_rotated;
3655
3656	//Buscar j i k que no canvin de r!!!
3657
3658	j=i-int(rho_pixel);
3659
3660	//cout<<"j inicial "<<j<<endl;
3661	int MINPIXinner[13]= {0,1,7,19,37,61,91,127,169,217,271,331,397};
3662	int MINPIXouter[5]={397,433,475,523,578};
3663
3664	if( r<12 ){
3665	if (j < MINPIXinner[r])
3666	{
3667	j = i+6 * r - int(rho_pixel);
3668	}
3669
3670	k=j-1;
3671	if (k < MINPIXinner[r])
3672	{
3673	k = MINPIXinner[r+1]-1;
3674	}
3675	}
3676	if( r > 11 && r < 16){
3677	if (j < MINPIXouter[r-12])
3678	{
3679	j = i + (r - 6) * 6 - int(rho_pixel);
3680
3681	}
3682
3683	k=j-1;
3684	if ( k < MINPIXouter[r-12])
3685	{
3686	k = MINPIXouter[r-11] - 1;
3687	}
3688	}
3689	rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k];
3690	}
3691	return (int(rho_pixel));
3692
3693	}
3694
3695
3696	//------------------------------------------------------------
3697	// @name produce_phes
3698	//
3699	// @desc read the photons of an event, pixelize them and simulate QE
3700	// @desc return various statistics and the array of Photoelectrons
3701	//
3702	// @date Mon Feb 14 16:44:21 CET 2000
3703	// @function @code
3704	//------------------------------------------------------------
3705
3706	int produce_phes( FILE *sp, // the input file
3707	class MGeomCam *camgeom, // the camera layout
3708	float minwl_nm, // the minimum accepted wavelength
3709	float maxwl_nm, // the maximum accepted wavelength
3710	class MTrigger *trigger, // the generated phes
3711	class MFadc *fadc,
3712	int *itotnphe, // total number of produced photoelectrons
3713	float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
3714	int *incph, // total number of cph read
3715	float *tmin_ns, // minimum arrival time of all phes
3716	float *tmax_ns, // maximum arrival time of all phes
3717	int ict // Telescope that is being analised to get the right QE.
3718	){
3719
3720	static uint i;
3721	static int k, ipixnum;
3722	static float cx, cy, wl, qe, t;
3723	static float cu, cv, cw;
3724	static MCCphoton photon;
3725	static float **qept;
3726	static char flag[SIZE_OF_FLAGS + 1];
3727	static float radius_mm;
3728
3729	// reset variables
3730
3731	for ( i=0; i<camgeom->GetNumPixels(); ++i ){
3732
3733	nphe[i] = 0.0;
3734
3735	}
3736
3737	*incph = 0;
3738
3739	radius_mm = camgeom->GetMaxRadius();
3740
3741
3742	//- - - - - - - - - - - - - - - - - - - - - - - - -
3743	// read photons and "map" them into the pixels
3744	//--------------------------------------------------
3745
3746	// initialize CPhoton
3747
3748	photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
3749
3750	// read the photons data
3751
3752
3753	// loop over the photons
3754
3755	while (1) {
3756
3757	fread ( flag, SIZE_OF_FLAGS, 1, sp );
3758
3759	if (isA( flag, FLAG_END_OF_EVENT ))
3760	break;
3761
3762	memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
3763
3764	fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
3765
3766
3767	// increase number of photons
3768
3769	(*incph)++;
3770
3771	// Chceck if photon is inside trigger time range
3772
3773	t = photon.get_t() ;
3774
3775	if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
3776	continue;
3777	}
3778
3779	//
3780	// Pixelization
3781	//
3782
3783	cx = photon.get_x();
3784	cy = photon.get_y();
3785
3786	// get wavelength
3787
3788	wl = photon.get_wl();
3789
3790	// cout << "wl " << wl << " radius " << sqrt(cxcx + cycy) << "\n";
3791
3792	// check if photon has valid wavelength and is inside outermost camera radius
3793
3794	if( (wl > maxwl_nm) \|\| (wl < minwl_nm) \|\| (sqrt(cxcx + cycy)*10 > radius_mm ) ){
3795	continue;
3796
3797	}
3798
3799	ipixnum = bpoint_is_in_pix(cx10, cy10, camgeom);
3800
3801	// -1 = the photon is in none of the pixels
3802	// 0 = the phton is in the central pixel, which is not used for trigger
3803	if (ipixnum==-1 \|\| ipixnum==0) {
3804	continue;
3805	}
3806
3807	//+++
3808	// QE simulation
3809	//---
3810
3811	// set pointer to the QE table of the relevant pixel
3812
3813	qept = (float **)QE[ict][ipixnum];
3814
3815	// check if wl is inside table; outside the table, QE is assumed to be zero
3816
3817	if((wl < qept[0][0]) \|\| (wl > qept[0][pointsQE[ict]-1])){
3818	continue;
3819
3820	}
3821
3822	// find data point in the QE table (-> k)
3823
3824	k = 1; // start at 1 because the condition was already tested for 0
3825	while (k < pointsQE[ict]-1 && qept[0][k] < wl){
3826	k++;
3827	}
3828
3829	// calculate the qe value between 0. and 1.
3830
3831	qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
3832
3833	// Apply incient angular correction due to Light Guides
3834
3835	cu=photon.get_u();
3836	cv=photon.get_v();
3837	cw=90.0-acos(sqrt(1-cucu-cvcv))*3.14159/180.0;
3838
3839	// find data point in the QE table (-> k)
3840
3841	k = 0;
3842	while (k < pointsWC-1 && WC[0][k] < cw){
3843	k++;
3844	}
3845
3846	// correct the qe with WC data.
3847
3848	qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
3849
3850	// if random > quantum efficiency, reject it
3851
3852	if ( (RandomNumber) > qe ) {
3853	continue;
3854	}
3855
3856	//+++
3857	// The photon has produced a photo electron
3858	//---
3859
3860	// cout << " accepted\n";
3861
3862	// increment the number of photoelectrons in the relevant pixel
3863
3864	nphe[ipixnum] += 1.0;
3865
3866	// store the new photoelectron
3867
3868	fadc->Fill(ipixnum,(t-*tmin_ns) ,
3869	trigger->FillShow(ipixnum,t-*tmin_ns),
3870	!((camgeom)[ipixnum].GetD()>(camgeom)[0].GetD()));
3871
3872	*itotnphe += 1;
3873	}
3874
3875	return(0);
3876
3877	}
3878
3879
3880	//------------------------------------------------------------
3881	// @name produce_nsbrates
3882	//
3883	// @desc read the starfield file, call produce_phes on it in,
3884	// @desc different wavebands, calculate the nsbrates
3885	//
3886	// @date Mon Feb 14 16:44:21 CET 2000
3887	// @function @code
3888	//------------------------------------------------------------
3889
3890	int produce_nsbrates( char *iname, // the starfield input file name
3891	MGeomCam *camgeom, // camera layout
3892	float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS], // the product of this function:
3893	// the NSB rates in phe/ns for each pixel
3894	int ict
3895	){
3896
3897	uint i, j; // counters
3898	int k, ii; // counters
3899
3900	static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
3901	WAVEBANDBOUND2,
3902	WAVEBANDBOUND3,
3903	WAVEBANDBOUND4,
3904	WAVEBANDBOUND5,
3905	WAVEBANDBOUND6 };
3906
3907	FILE *infile; // the input file
3908	fpos_t fileposition; // marker on the input file
3909	static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
3910	static MCEventHeader evth; // the event header
3911	static MCEventHeader evth_2; // the event header
3912	static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
3913	int reflector_file_version;
3914	int itnphe; // total number of produced photoelectrons
3915	int itotnphe; // total number of produced photoelectrons after averaging
3916	int incph; // total number of cph read
3917	float tmin_ns; // minimum arrival time of all phes
3918	float tmax_ns; // maximum arrival time of all phes
3919	float integtime_ns; // integration time from the starfield generator
3920	char flag_new[4];
3921
3922	// open input file
3923
3924	log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
3925
3926	infile = fopen( iname, "r" );
3927
3928	// check if the satrfield input file exists
3929
3930	if ( infile == NULL ) {
3931
3932	log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
3933
3934	log( SIGNATURE, "There is not NSB from the Stars\n");
3935
3936	return (0);
3937	}
3938
3939	// get signature, and check it
3940
3941	if((reflector_file_version=check_reflector_file( infile ))==FALSE){
3942	exit(1);
3943	}
3944
3945	// Instance of MTrigger and MFadc
3946
3947	MTrigger trigger(camgeom->GetNumPixels(),camgeom,Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
3948	MFadc flashadc(camgeom->GetNumPixels());
3949
3950	// initialize flag
3951
3952	strcpy( cflag, " \0" );
3953
3954	// get flag
3955
3956	fread( cflag, SIZE_OF_FLAGS, 1, infile );
3957
3958	if ( ! feof(infile)){
3959
3960	// reading .rfl file
3961
3962	if(!isA( cflag, FLAG_START_OF_RUN )){
3963	error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
3964	}
3965	else { // found start of run
3966
3967	if (reflector_file_version > 5){
3968
3969	fread( flag_new, 4, 1, infile );
3970
3971	if(!isA( flag_new, FLAG_START_OF_HEADER)){
3972
3973	// We break the main loop
3974	cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
3975	cout<<" We assume no Star Light"<<endl;
3976	return(0);
3977	}
3978
3979	Float_t flag_temp[SIZE_OF_MCRUNHEADER];
3980
3981	fread( flag_temp, (SIZE_OF_MCRUNHEADER)*sizeof(float), 1, infile );
3982
3983	}
3984
3985	fread( cflag, SIZE_OF_FLAGS, 1, infile );
3986
3987	if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
3988	fread( flag_new, 4, 1, infile );
3989
3990	if(!isA( flag_new, FLAG_EVENT_HEADER)){
3991
3992	// We break while events loop
3993	cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<endl;
3994	cout<<" We assume no light from Stars"<<endl;
3995	return(0);
3996	}
3997
3998	// get MCEventHeader
3999
4000	if (reflector_file_version<6)
4001	fread( (char*)&evth, evth.mysize(), 1, infile );
4002	else
4003	fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
4004
4005	if (reflector_file_version<6)
4006	integtime_ns = evth.get_energy();
4007	else
4008	integtime_ns = evth_2.get_energy();
4009
4010	// memorize where we are in the file
4011
4012	if (fgetpos( infile, &fileposition ) != 0){
4013	error( SIGNATURE, "Cannot position in file ...\n");
4014	}
4015
4016	// loop over the wavebands
4017
4018	for(i=0; i<iNUMWAVEBANDS; i++){
4019
4020	// initialize the rate array
4021
4022	for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4023	rate_phepns[j][i] = 0.;
4024	}
4025
4026	itotnphe = 0;
4027
4028	// read the photons and produce the photoelectrons
4029	// - in order to average over the QE simulation, call the
4030	// production function iNUMNSBPRODCALLS times
4031
4032	for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
4033
4034	// position the file pointer to the beginning of the photons
4035
4036	fsetpos( infile, &fileposition);
4037
4038	// produce photoelectrons
4039
4040	k = produce_phes( infile,
4041	camgeom,
4042	wl_nm[i],
4043	wl_nm[i+1],
4044	&trigger, // this is a dummy here
4045	&flashadc, // this is a dummy here
4046	&itnphe,
4047	nphe, // we want this!
4048	&incph,
4049	&tmin_ns,
4050	&tmax_ns,
4051	0
4052	); // photons from starfield
4053
4054	if( k != 0 ){ // non-zero returnvalue means error
4055	cout << "Exiting.\n";
4056	exit(1);
4057	}
4058
4059	for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4060	rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
4061	}
4062
4063	itotnphe += itnphe;
4064
4065	} // end for(ii=0 ...
4066
4067	fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
4068	wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
4069
4070	} // end for(i=0 ...
4071
4072	}
4073	else{
4074	cout << "Starfield file contains no event.\nExiting.\n";
4075	exit(1);
4076	} // end if( isA ... event
4077	} // end if ( isA ... run
4078	}
4079	else{
4080	cout << "Starfield file contains no run.\nExiting.\n";
4081	exit(1);
4082	}
4083
4084	fclose( infile );
4085
4086	return(0);
4087
4088	}
4089
4090
4091	//------------------------------------------------------------
4092	// @name produce_nsb_phes
4093	//
4094	// @desc produce the photoelectrons from the nsbrates
4095	//
4096	// @date Thu Feb 17 17:10:40 CET 2000
4097	// @function @code
4098	//------------------------------------------------------------
4099
4100	int produce_nsb_phes( float *atmin_ns,
4101	float *atmax_ns,
4102	float theta_rad,
4103	struct camera *cam,
4104	float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
4105	float difnsbr_phepns[iMAXNUMPIX],
4106	float extinction[iNUMWAVEBANDS],
4107	float fnpx[iMAXNUMPIX],
4108	MTrigger *trigger,
4109	MFadc *fadc,
4110	int *inphe,
4111	float base_mv[iMAXNUMPIX]){
4112
4113	float simtime_ns; // NSB simulation time
4114	int i, j, k, ii;
4115	float zenfactor; // correction factor calculated from the extinction
4116	int inumnsbphe; // number of photoelectrons caused by NSB
4117
4118	float t;
4119	TRandom random; // Random numbers generator
4120
4121	random.SetSeed(get_seeds(1));
4122
4123	ii = *inphe; // avoid dereferencing
4124
4125	// check if the arrival times are set; if not generate them
4126
4127	if(atmin_ns <SIMTIMEOFFSET_NS \|\| atmin_ns > *atmax_ns){
4128	*atmin_ns = 0.;
4129	*atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
4130
4131	}
4132	else{ // extend the event time window by the given offsets
4133
4134	atmin_ns = atmin_ns - SIMTIMEOFFSET_NS;
4135	atmax_ns = atmax_ns + SIMTIMEOFFSET_NS;
4136
4137	simtime_ns = atmax_ns - atmin_ns;
4138
4139	// make sure the simulated time is long enough for the FADC
4140	// simulation and not too long
4141
4142	if(simtime_ns< TOTAL_TRIGGER_TIME){
4143	atmin_ns = atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
4144	atmax_ns = atmin_ns + TOTAL_TRIGGER_TIME;
4145	simtime_ns = TOTAL_TRIGGER_TIME;
4146	}
4147
4148	if(simtime_ns> TOTAL_TRIGGER_TIME){
4149	atmax_ns =atmin_ns + TOTAL_TRIGGER_TIME;
4150	simtime_ns = TOTAL_TRIGGER_TIME;
4151	}
4152
4153	}
4154
4155	// initialize baselines
4156
4157	for(i=0; i<cam->inumpixels; i++){
4158	base_mv[i] = 0.;
4159	}
4160
4161	// calculate baselines and generate phes
4162
4163	for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
4164
4165	// calculate the effect of the atmospheric extinction
4166
4167	zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
4168
4169	for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
4170
4171	inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
4172	* simtime_ns );
4173
4174	base_mv[j] += inumnsbphe;
4175
4176	// randomize
4177
4178	if (inumnsbphe>0.0){
4179	inumnsbphe = random.Poisson (inumnsbphe );
4180	}
4181
4182	// create the photoelectrons
4183
4184	for(k=0; k < inumnsbphe; k++){
4185
4186	t=(RandomNumber * simtime_ns);
4187
4188	(fadc).Fill(j,t ,(trigger).FillNSB(j,t));
4189
4190	ii++; // increment total number of photoelectons
4191
4192	fnpx[j] += 1.; // increment number of photoelectrons in this pixel
4193
4194	}
4195
4196	} // end for(j=0 ...
4197	} // end for(i=0 ...
4198
4199	// finish baseline calculation
4200
4201	for(i=0; i<cam->inumpixels; i++){
4202	base_mv[i] = RESPONSE_FWHM RESPONSE_AMPLITUDE / simtime_ns;
4203	}
4204
4205	*inphe = ii; // update the pointer
4206
4207	return(0);
4208	}
4209
4210
4211	// @endcode
4212
4213
4214	//=------------------------------------------------------------
4215	//!@subsection Log of this file.
4216
4217	//!@{
4218	//
4219	// $Log: not supported by cvs2svn $
4220	// Revision 1.64 2003/10/21 07:42:50 blanch
4221	// A factor 2.35 to transform the fwhm into the sigma of gaussian was missing
4222	// in the storing of FADC single hpe pulse determination.
4223	//
4224	// Revision 1.63 2003/10/17 19:38:31 blanch
4225	// Now the camera program will stop if a undefined Geometry is required.
4226	// The NSB is internally scaled for any camera geometry and qe.
4227	// The seeds in the input card are used to initilize the random numbers.
4228	// The Amplitud stored in the MMcFadcHeader is the amplitud of the sphe reponse.
4229	// The Pedestal rms is simualted in a artificail empty event.
4230	//
4231	// Revision 1.62 2003/09/26 11:25:07 blanch
4232	// Modification to be able to read MGeomCam branch for any Geometry.
4233	//
4234	// Revision 1.61 2003/09/25 17:09:20 blanch
4235	// Bug on the number of phe from diffuse NSB fixed.
4236	//
4237	// Revision 1.60 2003/09/23 16:50:55 blanch
4238	// WE do not read ct_file anymore since all Telescope information is
4239	// in the reflector or in MGeomCam.
4240	//
4241	// Revision 1.58 2003/09/15 09:59:53 blanch
4242	// The concept of the camera prgoram has not changed but this version has
4243	// quite a lot of changes to allow several Camera geometries as well as
4244	// multitelescope.
4245	//
4246	// It is suposed to be the first working code for camera 0.7.
4247	//
4248	// Revision 1.57 2003/07/17 18:02:46 blanch
4249	// Several new features introduced as well as fixed bugs
4250	//
4251	// - 1/100 events printed out
4252	// - Low gain implemented
4253	// - Different response for outer and inner pixels
4254	// - Some warnings removed
4255	// - pedestal and qe file from inpuit card
4256	// - Faster electronic simulation
4257	//
4258	// Revision 1.55 2003/02/12 12:22:10 blanch
4259	// * empty log message *
4260	//
4261	// Revision 1.54 2003/02/12 11:55:01 blanch
4262	// * empty log message *
4263	//
4264	// Revision 1.53 2003/01/23 18:35:21 blanch
4265	// * empty log message *
4266	//
4267	// Revision 1.52 2003/01/20 17:19:20 blanch
4268	// It fills the MMcCorsikaRun.
4269	//
4270	// Revision 1.51 2003/01/14 13:40:17 blanch
4271	// MRawRunHeader::fNumEvents has been filled with the number of events in
4272	// this file.
4273	// Problems in fImpact computation have been solved.
4274	// Option to set a dc value to rise the discriminator threshold has been added.
4275	//
4276	// Revision 1.50 2003/01/07 16:33:31 blanch
4277	// Star Field Rotation has been implemented by Raul Orduna. Now there is a
4278	// rotation for each shower. It is done by a non enter pixel rotation assuming
4279	// a circular symetry of the camera. It is not exact but it is accurate enough and
4280	// much faster.
4281	//
4282	// Revision 1.49 2002/12/13 10:04:07 blanch
4283	// * empty log message *
4284	//
4285	// Revision 1.48 2002/12/12 17:40:50 blanch
4286	// * empty log message *
4287	//
4288	// Revision 1.47 2002/12/10 17:19:31 blanch
4289	// * empty log message *
4290	//
4291	// Revision 1.46 2002/11/08 17:51:00 blanch
4292	// * empty log message *
4293	//
4294	// Revision 1.45 2002/10/29 17:15:27 blanch
4295	// Reading several reflector versions.
4296	//
4297	// Revision 1.44 2002/10/18 16:53:03 blanch
4298	// Modification to read several reflector files.
4299	//
4300	// Revision 1.43 2002/09/13 10:53:39 blanch
4301	// Minor change to remove some undisired comments.
4302	//
4303	// Revision 1.42 2002/09/09 16:00:49 blanch
4304	// Statement has been included to avoid writting to disk MParContainer and MArray.
4305	// It has also been added the effect of the WC, the actual values must be added,
4306	// once they are measured.
4307	//
4308	// Revision 1.41 2002/09/04 09:57:42 blanch
4309	// Modifications done to use MGeomCam from MARS.
4310	//
4311	// Revision 1.40 2002/07/16 16:15:22 blanch
4312	// A first implementation of the Star field rotation has been done.
4313	//
4314	// Revision 1.39 2002/04/27 10:48:39 blanch
4315	// Some unused varibles have been removed.
4316	//
4317	// Revision 1.38 2002/03/18 18:44:29 blanch
4318	// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
4319	//
4320	// Revision 1.37 2002/03/18 16:42:20 blanch
4321	// The data member fProductionSite of the MMcRunHeader has been set to 0,
4322	// which means that the prodution site is unknown.
4323	//
4324	// Revision 1.35 2002/03/15 15:15:52 blanch
4325	// Several modifications to simulate the actual trigger zone.
4326	//
4327	// Revision 1.34 2002/03/13 18:13:56 blanch
4328	// Some changes to fill correctly the new format of MMcRunHeader.
4329	//
4330	// Revision 1.33 2002/03/04 17:21:48 blanch
4331	// Small and not important changes.
4332	//
4333	// Revision 1.32 2002/02/28 15:04:52 blanch
4334	// A small back has been solved. Before, while not using the option
4335	// writte_all_images, not all triggered showers were stored. Now it is solved.
4336	// For that it is importatn taht the less restrictive trigger option is
4337	// checked first.
4338	// A new facility has been introduced and now one can choose the step size in
4339	// trigger loop mode for the discriminator threshold.
4340	// The close-compact topology for less than 3 pixels does not make sense. Before
4341	// the program was ignoring that, now it switch to simple neighbour condition.
4342	//
4343	// Revision 1.31 2002/01/18 17:41:02 blanch
4344	// The option of adding noise to all pixels or to not adding the noise
4345	// has been added.
4346	// We removed the pixels larger than 577. When there were more than one
4347	// trigger in one shower, the pixel number was increasing. Now it is
4348	// flagged by the variable MMcTrig::fFirstLvlTrig.
4349	//
4350	// Revision 1.30 2001/11/27 09:49:54 blanch
4351	// Fixing bug which was treating wrongly the extension of star photons.
4352	//
4353	// Revision 1.29 2001/11/14 17:38:23 blanch
4354	// Sveral changes have been done:
4355	// - bpoint_in_in_pixel has been dodified to speed up the program
4356	// - Some minor changes have been done to remove warnings
4357	// - buffer size and split version of the Branches have been removed
4358	// - Some modifications were needed to adat the program to the new
4359	// MRawEvtData::DeletePixels
4360	//
4361	// Revision 1.28 2001/10/26 16:31:45 blanch
4362	// Removing several warnings.
4363	//
4364	// Revision 1.27 2001/09/05 10:04:33 blanch
4365	// * empty log message *
4366	//
4367	// Revision 1.26 2001/07/19 09:29:53 blanch
4368	// Different threshold for each pixel can be used.
4369	//
4370	// Revision 1.25 2001/05/08 08:07:54 blanch
4371	// New numbering for branches from different trigger conditions has been
4372	// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
4373	// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
4374	// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
4375	// information already filled.
4376	//
4377	// Revision 1.24 2001/04/06 16:48:09 magicsol
4378	// New camera version able to read the new format of the reflector output:
4379	// reflector 0.4
4380	//
4381	// Revision 1.23 2001/03/28 16:13:41 blanch
4382	// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
4383	// been solved.
4384	//
4385	// Revision 1.22 2001/03/20 18:52:43 blanch
4386	// * empty log message *
4387	//
4388	// Revision 1.21 2001/03/19 19:53:03 blanch
4389	// Some print outs have been removed.
4390	//
4391	// Revision 1.20 2001/03/19 19:30:06 magicsol
4392	// Minor changes have been done to improve the FADC pedestals treatment.
4393	//
4394	// Revision 1.19 2001/03/05 11:14:41 magicsol
4395	// I changed the position of readinf a parameter. It is a minnor change.
4396	//
4397	// Revision 1.18 2001/03/05 10:36:52 blanch
4398	// A branch with information about the FADC simulation (MMcFadcHeader) is writen
4399	// in the McHeaders tree of the aoutput root file.
4400	// The NSB contribution is only applied if the the number of phe form the shower
4401	// are above a given number.
4402	//
4403	// Revision 1.17 2001/02/23 11:05:57 magicsol
4404	// Small changes due to slightly different output format and the introduction of
4405	// pedesals for teh FADC.
4406	//
4407	// Revision 1.16 2001/01/18 18:44:40 magicsol
4408	// * empty log message *
4409	//
4410	// Revision 1.15 2001/01/17 09:32:27 harald
4411	// The changes are neccessary to have the same name for trees in MC and in
4412	// data. So now there should be no differences in MC data and real data from
4413	// FADC system.
4414	//
4415	// Revision 1.14 2001/01/15 12:33:34 magicsol
4416	// Some modifications have been done to use the new (Dec'2000) Raw data format.
4417	// There are also some minnors modifications to adapt the improvements in the
4418	// MTrigger class (overlaping time and trigger cells).
4419	//
4420	// Revision 1.13 2000/10/25 08:14:23 magicsol
4421	// The routine that produce poisson random numbers to decide how many phe
4422	// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
4423	// is used.
4424	//
4425	// Revision 1.10 2000/07/04 14:10:20 MagicSol
4426	// Some changes have been done in the root output file. The RawEvt tree is only
4427	// stored in the single trigger mode.
4428	// The trigger input parameters are also given by the general input card.
4429	// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
4430	// each one can be studied seperately.
4431	//
4432	// Revision 1.9 2000/06/13 13:25:24 blanch
4433	// The multiple arrays have been replaced, since they do not work
4434	// in alpha machines. Now we are using pointers and the command new
4435	// to allocate memory.
4436	//
4437	// Revision 1.8 2000/05/11 13:57:27 blanch
4438	// The option to loop over several trigger configurations has been included.
4439	// Some non-sense with triggertime range has been solved.
4440	// Montecarlo information and ADC counts are saved in a root file.
4441	// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
4442	//
4443	// Revision 1.7 2000/03/24 18:10:46 blanch
4444	// A first FADC simulation and a trigger simulation are already implemented.
4445	// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
4446	// 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.
4447	//
4448	// Revision 1.6 2000/03/20 18:35:11 blanch
4449	// 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.
4450	//
4451	// Revision 1.5 2000/02/18 17:40:35 petry
4452	// This version includes drastic changes compared to camera.cxx 1.4.
4453	// It is not yet finished and not immediately useful because the
4454	// trigger simulation is not yet re-implemented. I had to take it
4455	// out together with some other stuff in order to tidy the whole
4456	// program up. This is not meant as an insult to anyone. I needed
4457	// to do this in order to be able to work on it.
4458	//
4459	// This version has been put in the repository in order to be
4460	// able to share the further development with others.
4461	//
4462	// If you need something working, wait or take an earlier one.
4463	// See file README.
4464	//
4465	// Revision 1.4 2000/01/25 08:36:23 petry
4466	// The pixelization in previous versions was buggy.
4467	// This is the first version with a correct pixelization.
4468	//
4469	// Revision 1.3 2000/01/20 18:22:17 petry
4470	// Found little bug which makes camera crash if it finds a photon
4471	// of invalid wavelength. This bug is now fixed and the range
4472	// of valid wavelengths extended to 290 - 800 nm.
4473	// This is in preparation for the NSB simulation to come.
4474	// Dirk
4475	//
4476	// Revision 1.2 1999/11/19 08:40:42 harald
4477	// Now it is possible to compile the camera programm under osf1.
4478	//
4479	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
4480	// This the starting point for CVS controlled further developments of the
4481	// camera program. The program was originally written by Jose Carlos.
4482	// But here you can find a "rootified" version to the program. This means
4483	// that there is no hbook stuff in it now. Also the output of the
4484	// program changed to the MagicRawDataFormat.
4485	//
4486	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
4487	//
4488	// Revision 1.3 1999/10/22 15:01:28 petry
4489	// version sent to H.K. and N.M. on Fri Oct 22 1999
4490	//
4491	// Revision 1.2 1999/10/22 09:44:23 petry
4492	// first synthesized version which compiles and runs without crashing;
4493	//
4494	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
4495	// first synthesised version
4496	//
4497	// Revision 1.13 1999/03/15 14:59:05 gonzalez
4498	// camera-1_1
4499	//
4500	// Revision 1.12 1999/03/02 09:56:10 gonzalez
4501	// * empty log message *
4502	//
4503	//
4504	//!@}
4505
4506	//=EOF

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

Download in other formats: