Context Navigation

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

Visit:

Last change on this file since 2796 was 2453, checked in by blanch, 22 years ago
Independent reference time for each telescope.
File size: 128.2 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.66 $
24	// $Author: blanch $
25	// $Date: 2003-11-03 12:58:40 $
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.SetSeed(Int_t(get_seeds(1)/get_seeds(0))+1);
1322
1323	lons.SetPath(nsbpathname);
1324
1325	// Instance of the Mlons class
1326	MLons lons_outer(Trigger_response_ampl, Trigger_response_fwhm,
1327	FADC_resp_ampl_out,FADC_resp_fwhm_out);
1328
1329	lons_outer.SetSeed(Int_t(get_seeds(1)/get_seeds(0))+2);
1330
1331	lons_outer.SetPath(nsbpath_outer);
1332
1333	if( simulateNSB ){
1334
1335	//
1336	// Calculate the non-diffuse NSB photoelectron rates
1337	//
1338
1339	// FIXME --- star NSB different for each camera?
1340	log(SIGNATURE,"Produce NSB rates from Star Field");
1341
1342	k = produce_nsbrates( starfieldname,
1343	((MGeomCam*)(camgeom.UncheckedAt(0))),
1344	nsbrate_phepns,
1345	0);
1346
1347	if (k != 0){
1348	cout << "Error when reading starfield... \nExiting.\n";
1349	exit(1);
1350	}
1351
1352	// calculate diffuse rate correcting for the pixel size and telescope
1353
1354	float factorNSB_ct;
1355	for(int ict=0;ict<ct_Number;ict++){
1356
1357	// First we set the factor due to the mirror size respect to the MAGIC mirror.
1358
1359	switch(GeometryCamera[ict]){
1360	case 1:
1361	case 2:
1362	case 3:
1363	case 4:
1364	factorNSB_ct=1.0;
1365	break;
1366	case 5:
1367	case 6:
1368	case 7:
1369	factorNSB_ct=1000.0/239.0;
1370	break;
1371	case 8:
1372	case 9:
1373	default:
1374	factorNSB_ct=1.0;
1375	break;
1376	}
1377
1378	factorNSB_ct=factorNSB_ct/
1379	((((MGeomCam)(camgeom.UncheckedAt(ict)))).GetCameraDist()1000
1380	(((MGeomCam)(camgeom.UncheckedAt(ict)))).GetCameraDist()1000
1381	PIXEL_SIZE*PIXEL_SIZE);
1382
1383	for(UInt_t ui=0;
1384	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); ui++){
1385	const Float_t size=
1386	(((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD();
1387	diffnsb_phepns[ict][ui] = (Int_t(meanNSBfactorNSB_ct100sizesize+0.5))/(100.0)*factorqe_NSB[ict];
1388	}
1389	}
1390
1391	// calculate nsb rate including diffuse and starlight
1392	// we also include the extinction effect
1393	for(int ict=0;ict<ct_Number;ict++){
1394	for(j=0;j<iNUMWAVEBANDS;j++){
1395	// calculate the effect of the atmospheric extinction
1396
1397	zenfactor = pow(10., -0.4 * ext[j] );
1398
1399	for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1400	nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
1401	nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
1402	}
1403	}
1404	}
1405	}
1406
1407
1408	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1409	//
1410	// Now 500 empty event with the conditin in which the camera is run
1411	// are simulated. IN this way one gets an estimation of the
1412	// sigma for the pedestal of each FADC channel.
1413	// This sigma computtaion is done assuming any noise taht effects
1414	// the FADC but there is no rotation of the Star Field (otherwise it
1415	// should be done for each event)
1416	//
1417	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1418
1419	for(int ict=0;ict<ct_Number;ict++){
1420	for(UInt_t ui=0;
1421	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1422	ui++){
1423	fadc_sigma[ui]=0.0;
1424	fadc_sigma_low[ui]=0.0;
1425	}
1426	for(int ie=0;ie<500;ie++){
1427	Fadc_CT[ict]->Reset() ;
1428	if (addElecNoise){
1429	Fadc_CT[ict]->ElecNoise(FADC_noise) ;
1430	}
1431	if(simulateNSB){
1432	for(UInt_t ui=0;
1433	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1434	ui++){
1435	if(nsb_phepns[ict][ui]>0.0){
1436	if((((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD()>(((MGeomCam)(camgeom.UncheckedAt(ict))))[0].GetD()){
1437	k=lons_outer.GetResponse(nsb_phepns[ict][ui],0.01,
1438	& nsb_trigresp[0],
1439	& nsb_fadcresp[0]);
1440	}
1441	else{
1442	k=lons.GetResponse(nsb_phepns[ict][ui],0.01,
1443	& nsb_trigresp[0],& nsb_fadcresp[0]);
1444	}
1445	if(k==0){
1446	cout << "Exiting.\n";
1447	exit(1);
1448	}
1449	Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
1450	}
1451	}
1452	}
1453	Fadc_CT[ict]->Pedestals();
1454	for(UInt_t ui=0;
1455	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1456	ui++){
1457	fadc_sigma[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,1))/500.0;
1458	fadc_sigma_low[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,0))/500.0;
1459	}
1460	}
1461	HeaderFadc[ict]->SetPedestalSigma(&fadc_sigma_low[0],&fadc_sigma[0],
1462	((MGeomCam*)(camgeom.UncheckedAt(ict)))
1463	->GetNumPixels());
1464	}
1465
1466	//
1467	// Read the reflector file with the Cherenkov data
1468	//
1469
1470	// select input file
1471
1472	if ( Data_From_STDIN ) {
1473
1474	inputfile[0] = stdin;
1475
1476	}
1477	else{
1478
1479	for(int ict=0;ict<ct_Number;ict++){
1480
1481	log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
1482	inputfile[ict] = fopen( inname_CT[ict], "r" );
1483
1484	if ( inputfile[ict] == NULL )
1485	error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
1486	}
1487
1488	}
1489
1490	// get signature, and check it
1491
1492	for(int ict=0;ict<ct_Number;ict++){
1493	if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE){
1494	exit(1);
1495	}
1496	}
1497
1498	// open data file
1499
1500	log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
1501	datafile.open( datname );
1502
1503	if ( datafile.bad() )
1504	error( SIGNATURE, "Cannot open data file: %s\n", datname );
1505
1506	// initializes flag
1507
1508	strcpy( flag, " \0" );
1509	strcpy( flag_new, " \0" );
1510
1511	// allocate space for PMTs numbers of pixels
1512
1513	fnpix = new float [ct_NPixels];
1514
1515
1516	//!@' @#### Main loop.
1517	//@'
1518
1519	// get flag
1520
1521	for(int ict=0;ict<ct_Number;ict++)
1522	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1523
1524
1525	// loop over the file
1526
1527	still_in_loop = TRUE;
1528
1529	// FIXME --- check if not eof for all input reflector files
1530
1531	while (
1532	((! Data_From_STDIN) && ( !feof(inputfile[0]) ))
1533	\|\|
1534	(Data_From_STDIN && still_in_loop)
1535	) {
1536
1537	// reading .rfl files
1538	if(!isA( flag, FLAG_START_OF_RUN )){
1539
1540	// We break the main loop
1541	cout<<"Warning: Expected start of run flag, but found:"<<flag<<endl;
1542	cout<<" We break the main loop"<<endl;
1543	break;
1544	}
1545	else { // found start of run
1546
1547	for(int ict=0;ict<ct_Number;ict++) {
1548
1549	fread( flag_new, 4, 1, inputfile[ict] );
1550
1551	if(!isA( flag_new, FLAG_START_OF_HEADER)){
1552
1553	// We break the main loop
1554	cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
1555	cout<<" We break the main loop"<<endl;
1556	break;
1557	}
1558
1559	fread((char)&mcrunh,(SIZE_OF_MCRUNHEADER)sizeof(float), 1, inputfile[ict] );
1560	}
1561
1562	nshow=0;
1563
1564	for(int ict=0;ict<ct_Number;ict++){
1565	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1566	}
1567	while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
1568	for(int ict=0;ict<ct_Number;ict++) {
1569	fread( flag_new, 4, 1, inputfile[ict] );
1570	if(!isA( flag_new, FLAG_EVENT_HEADER)){
1571
1572	// We break while events loop
1573	cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<" "<<ict<<endl;
1574	cout<<" We break the while events loop"<<endl;
1575	break;
1576	}
1577	}
1578
1579
1580	//
1581	// Clear Trigger and Fadc
1582	//
1583	for(int ict=0;ict<ct_Number;ict++) {
1584	Trigger_CT[ict]->Reset() ;
1585	Trigger_CT[ict]->ClearFirst();
1586	Trigger_CT[ict]->ClearZero();
1587	Fadc_CT[ict]->Reset() ;
1588	}
1589
1590	++nshow;
1591	if((nshow+ntshow)%100 == 1)
1592	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
1593
1594	// get MCEventHeader
1595	if (reflector_file_version<6){
1596	for(int ict=0;ict<ct_Number;ict++)
1597	fread( (char*)&mcevth[ict], mcevth[ict].mysize(),
1598	1, inputfile[ict] );
1599	}
1600	else{
1601	for(int ict=0;ict<ct_Number;ict++)
1602	fread( (char*)&mcevth_2[ict], mcevth_2[ict].mysize(),
1603	1, inputfile[ict] );
1604	}
1605
1606	// calculate impact parameter (shortest distance betwee the original
1607	// trajectory of the primary (assumed shower-axis) and the
1608	// direction where the telescope points to
1609	//
1610	// we use the following equation, given that the shower core position
1611	// is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
1612	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
1613	// and (l2,m2,n2)
1614	//
1615	// \| \|
1616	// \| x1-x2 y1-y2 z1-z2 \|
1617	// \| \|
1618	// + \| l1 m1 n1 \|
1619	// - \| \|
1620	// \| l2 m2 n2 \|
1621	// \| \|
1622	// dist = ------------------------------------ ( > 0 )
1623	// [ \|l1 m1\|2 \|m1 n1\|2 \|n1 l1\|2 ]1/2
1624	// [ \| \| + \| \| + \| \| ]
1625	// [ \|l2 m2\| \|m2 n2\| \|n2 l2\| ]
1626	//
1627	// playing a little bit, we get this reduced for in our case:
1628	//
1629	//
1630	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
1631	// [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
1632	// 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
1633
1634	// read the direction of the incoming shower
1635	// It is done only for one telescope since it is suposed
1636	// to be the same shower for all of them
1637	if (reflector_file_version<6){
1638	thetashw = mcevth[0].get_theta();
1639	phishw = mcevth[0].get_phi();
1640	}
1641	else{
1642	thetashw = mcevth_2[0].get_theta();
1643	phishw = mcevth_2[0].get_phi();
1644	}
1645
1646	// calculate vector for shower
1647
1648	l1 = sin(thetashw)*cos(phishw);
1649	m1 = sin(thetashw)*sin(phishw);
1650	n1 = cos(thetashw);
1651
1652	if (reflector_file_version<6){
1653
1654	mcevth[0].get_core(&coreX, &coreY);
1655
1656	for(int ict=0;ict<ct_Number;ict++){
1657	// read the deviation of the telescope with respect to the shower
1658	mcevth[ict].get_deviations ( &thetaCT[ict], &phiCT[ict] );
1659
1660	if ( (thetaCT[ict] == 0.) && (phiCT[ict] == 0.) ) {
1661
1662	// CT was looking to the source (both lines are parallel)
1663	// therefore, we calculate the impact parameter as the distance
1664	// between the CT axis and the core position
1665
1666	impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1667	thetaCT[ict] += thetashw;
1668	phiCT[ict] += phishw;
1669	} else {
1670
1671	// the shower comes off-axis
1672
1673	// obtain with this the final direction of the CT
1674
1675	thetaCT[ict] += thetashw;
1676	phiCT[ict] += phishw;
1677
1678	// calculate vector for telescope
1679
1680	l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
1681	m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
1682	n2 = cos(thetaCT[ict]);
1683
1684	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1685	den = (SQR(l1m2 - l2m1) +
1686	SQR(m1n2 - m2n1) +
1687	SQR(n1l2 - n2l1));
1688	den = sqrt(den);
1689	impactD[ict] = fabs(num)/den;
1690	}
1691	}
1692	}
1693	else{
1694	mcevth_2[0].get_core(&coreX, &coreY);
1695	for(int ict=0;ict<ct_Number;ict++){
1696	thetaCT[ict]=mcevth_2[ict].get_theta_CT();
1697	phiCT[ict]=mcevth_2[ict].get_phi_CT();
1698	if ( (thetaCT[ict] == thetashw) && (phiCT[ict] == phishw) ) {
1699
1700	// CT was looking to the source (both lines are parallel)
1701	// therefore, we calculate the impact parameter as the distance
1702	// between the CT axis and the core position
1703
1704	impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1705
1706	} else {
1707
1708	// the shower comes off-axis
1709
1710	// calculate vector for telescope
1711
1712	l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
1713	m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
1714	n2 = cos(thetaCT[ict]);
1715
1716	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1717	den = (SQR(l1m2 - l2m1) +
1718	SQR(m1n2 - m2n1) +
1719	SQR(n1l2 - n2l1));
1720	den = sqrt(den);
1721
1722	impactD[ict] = fabs(num)/den;
1723
1724	}
1725	}
1726	}
1727
1728	// energy cut
1729
1730	if ( Select_Energy ) {
1731	if (reflector_file_version<6)
1732	if (( mcevth[0].get_energy() < Select_Energy_le ) \|\|
1733	( mcevth[0].get_energy() > Select_Energy_ue )) {
1734	log(SIGNATURE, "select_energy: shower rejected.\n");
1735	continue;
1736	}
1737	else
1738	if (( mcevth_2[0].get_energy() < Select_Energy_le ) \|\|
1739	( mcevth_2[0].get_energy() > Select_Energy_ue )) {
1740	log(SIGNATURE, "select_energy: shower rejected.\n");
1741	continue;
1742	}
1743	}
1744
1745	inumphe=0;
1746
1747	for(int ict=0;ict<ct_Number;ict++){
1748
1749	// Read first and last time and put inumphe_CT[0] to 0
1750
1751	if (reflector_file_version<6)
1752	mcevth[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1753	else
1754	mcevth_2[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1755
1756	inumphe_CT[ict]=0;
1757
1758	// read the photons and produce the photoelectrons
1759
1760	k = produce_phes( inputfile[ict],
1761	((MGeomCam*)(camgeom.UncheckedAt(ict))),
1762	WAVEBANDBOUND1,
1763	WAVEBANDBOUND6,
1764	Trigger_CT[ict], // will be changed by the function!
1765	Fadc_CT[ict], // will be changed by the function!
1766	&inumphe_CT[ict], // important for later: the size of photoe[]
1767	fnpix, // will be changed by the function!
1768	&ncph[ict], // will be changed by the function!
1769	&arrtmin_ns, // will be changed by the function!
1770	&arrtmax_ns, // will be changed by the function!
1771	ict);
1772
1773	inumphe=(inumphe<inumphe_CT[ict])?inumphe_CT[ict]:inumphe;
1774
1775	if( k != 0 ){ // non-zero returnvalue means error
1776	cout << "Exiting.\n";
1777	exit(1);
1778	}
1779	}
1780
1781	// NSB simulation
1782
1783	if(simulateNSB && nphe2NSB<=inumphe){
1784
1785	if(Starfield_rotate){
1786
1787	// Introduction rho angle
1788
1789	zenith = thetashw;
1790	azimutal = phishw;
1791	C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
1792	C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
1793	C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
1794	rho = acos( C1/C2 );
1795
1796	if ( sin(azimutal) < 0)
1797	{
1798	rho = 2 * 3.14159 - rho;
1799	}
1800	else
1801	{
1802	rho = rho;
1803	}
1804	rho = rho*180/3.14159;
1805
1806	// Rottion of the NSB
1807	// FIXME --- We should rotate for all cameras. Is it always the same rho?
1808	for(int ict=0;ict<ct_Number;ict++){
1809	k = size_rotated(
1810	&nsb_phepns_rotated[ict][0],
1811	nsb_phepns[ict],
1812	rho);
1813	}
1814	}
1815
1816	// Fill trigger and fadc response in the trigger class from the database
1817	for(int ict=0;ict<ct_Number;ict++){
1818
1819	for(UInt_t ui=0;ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1820	if(nsb_phepns_rotated[ict][ui]>0.0){
1821	if((((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD()>(((MGeomCam)(camgeom.UncheckedAt(ict))))[0].GetD()){
1822	k=lons_outer.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
1823	& nsb_trigresp[0],
1824	& nsb_fadcresp[0]);
1825	}
1826	else{
1827	k=lons.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
1828	& nsb_trigresp[0],& nsb_fadcresp[0]);
1829	}
1830	if(k==0){
1831	cout << "Exiting.\n";
1832	exit(1);
1833	}
1834	Trigger_CT[ict]->AddNSB(ui,nsb_trigresp);
1835	Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
1836	}
1837	}
1838	}
1839
1840	}// end if(simulateNSB && nphe2NSB<=inumphe_CT[0]) ...
1841
1842	for(int ict=0;ict<ct_Number;ict++){
1843	inumphensb[ict]=0;
1844	for (UInt_t ui=0;ui<
1845	((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1846	inumphensb[ict]+=nsb_phepns[ict][ui]*TOTAL_TRIGGER_TIME;
1847	}
1848	ntcph[ict]+=ncph[ict];
1849	if ((nshow+ntshow)%100 == 1){
1850	log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1851	ncph[ict], ntcph[ict]);
1852
1853	cout << "Total number of phes in CT "<<ict<<": "
1854	<< inumphe_CT[ict]<<" (+ ";
1855	cout<<inumphensb[ict]<<" mean expected number from NSB)"<<endl;
1856	}
1857	}
1858
1859	// skip it ?
1860
1861	int i;
1862	for (i=0; i<nSkip; ++i ) {
1863	if (Skip[i] == (nshow+ntshow)) {
1864	i = -1;
1865	break;
1866	}
1867	}
1868
1869	// if after the previous loop, the exit value of i is -1
1870	// then the shower number is in the list of showers to be
1871	// skipped
1872
1873	if (i == -1) {
1874	log(SIGNATURE, "\t\tskipped!\n");
1875	continue;
1876	}
1877
1878	//++++++++++++++++++++++++++++++++++++++++++++++++++
1879	// at this point we have a camera full of
1880	// ph.e.s
1881	// we should first apply the trigger condition,
1882	// and if there's trigger, then clean the image,
1883	// calculate the islands statistics and the
1884	// other parameters of the image (Hillas' parameters
1885	// and so on).
1886	//--------------------------------------------------
1887
1888	// TRIGGER HERE
1889
1890	// We should simulate the AC coupling behaviour:
1891	// For the FADC it is only done for the NSB while producinr
1892	// the StarResponse database.
1893	// For the trigger is done in the Trigger.Diskriminate(), which
1894	// is called later (it should be seperated to speed up the program
1895	//
1896
1897	// now the noise of the electronic
1898	// (preamps, optical transmission,..) is introduced.
1899	// This is done inside the class MTrigger by the method ElecNoise.
1900	//
1901
1902	for(int ict=0;ict<ct_Number;ict++) {
1903	if (addElecNoise && nphe2NSB<=inumphe){
1904
1905	Trigger_CT[ict]->ElecNoise(Trigger_noise) ;
1906
1907	Fadc_CT[ict]->ElecNoise(FADC_noise) ;
1908	}
1909	}
1910
1911	// now a shift in the fadc signal due to the pedestlas is
1912	// intorduced
1913	// This is done inside the class MFadc by the method Pedestals
1914	for(int ict=0;ict<ct_Number;ict++) {
1915	Fadc_CT[ict]->Pedestals();
1916	}
1917
1918	// We study several trigger conditons
1919	if(Trigger_Loop){
1920
1921	// Set to zero the flag to know if some conditon has triggered
1922	btrigger=0;
1923	flagstoring = 0;
1924
1925	// Loop over trigger threshold
1926	int iconcount;
1927	for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1928	for (int i=0;i<ct_NPixels;i++){
1929	fpixelthres[i]=
1930	((Float_t)(fthrescount)>=qThreshold[0][i])?
1931	(Float_t)(fthrescount):qThreshold[0][i];
1932
1933	// Rise the discrimnator threshold to avoid huge rates
1934
1935	if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
1936	for(int ii=0;ii<ct_NPixels;ii++)
1937	if( nsb_phepns_rotated[0][ii]>riseDiskThres)
1938	fpixelthres[ii]=secureDiskThres;
1939	}
1940	Trigger_CT[0]->SetThreshold(fpixelthres);
1941
1942	Trigger_CT[0]->Diskriminate();
1943
1944
1945	//
1946	// look if in all the signals in the trigger signal branch
1947	// is a possible Trigger. Therefore we habe to diskriminate all
1948	// the simulated analog signals (Method Diskriminate in class
1949	// MTrigger). We look simultanously for the moments at which
1950	// there are more than TRIGGER_MULTI pixels above the
1951	// CHANNEL_THRESHOLD.
1952	//
1953
1954	// Set trigger flags to zero
1955	Lev0=0;
1956	Lev1=0;
1957
1958	// loop over multiplicity of trigger configuration
1959	for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1960	Trigger_CT[0]->SetMultiplicity(imulticount);
1961	Trigger_CT[0]->ClearZero();
1962
1963	Lev0=(Short_t) Trigger_CT[0]->ZeroLevel();
1964	if (Lev0>0 \|\| Write_All_Images \|\| btrigger){
1965
1966	// loop over topologies
1967	for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1968	Lev1=0;
1969
1970	if(itopocount==0 && imulticount>7) continue;
1971	//COBB if(itopocount==2 && imulticount<3) continue;
1972	// It only makes to look for a different topology
1973	// if there are 3 or more N pixels.
1974	if(imulticount<3)
1975	Trigger_CT[0]->SetTopology(1);
1976	else
1977	{
1978	// We should be careful that topologies are sort from
1979	// the less to the more restrictive one.
1980	Trigger_CT[0]->SetTopology(isorttopo[itopocount]);
1981	}
1982	Trigger_CT[0]->ClearFirst();
1983
1984	//
1985	// Start the First Level Trigger simulation
1986	//
1987	if(Lev0!=0)
1988	Lev1=Trigger_CT[0]->FirstLevel();
1989	if(Lev1>0) {
1990	btrigger= 1;
1991	ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++;
1992	}
1993
1994	Lev0=1;
1995	Int_t NumImages = Lev1;
1996	if(Lev1==0 && (Write_All_Images \|\| btrigger)){
1997	btrigger= 1;
1998	NumImages=1;
1999	Lev0=0;
2000	}
2001
2002	for (Int_t ii=0;ii<NumImages;ii++){
2003	if (Write_McTrig){
2004	McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
2005	McTrig[iconcount]->SetTime(Trigger_CT[0]->GetFirstLevelTime(ii),ii+1);
2006	Trigger_CT[0]->GetMapDiskriminator(trigger_map);
2007	McTrig[iconcount]->SetMapPixels(trigger_map,ii);
2008	}
2009	//
2010	// fill inside the class fadc the member output
2011	//
2012
2013	Fadc_CT[0]->TriggeredFadc(Trigger_CT[0]->GetFirstLevelTime(ii));
2014
2015	if( Write_RawEvt ){
2016	//
2017	// Fill the header of this event
2018	//
2019
2020	EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
2021	// fill pixel information
2022	if (Lev1 \|\| Write_All_Images){
2023	if (addElecNoise) Fadc_CT[0]->DigitalNoise();
2024	for(UInt_t i=0;
2025	i<((MGeomCam*)(camgeom.UncheckedAt(0)))
2026	->GetNumPixels();i++){
2027	if(!Fadc_CT[0]->IsPixelUsed(i)) continue;
2028	for (j=0;j<FADC_SLICES;j++){
2029	fadcValues->AddAt(Fadc_CT[0]->GetFadcSignal(i,j),j);
2030	fadcValuesLow->AddAt(Fadc_CT[0]->GetFadcLowGainSignal(i,j),j);
2031	}
2032	EvtData[iconcount]->AddPixel(i,fadcValues,0);
2033	EvtData[iconcount]->AddPixel(i,fadcValuesLow,kTRUE);
2034	}
2035	}
2036	}
2037	}
2038	//
2039	// Increase counter of analised trigger conditions
2040	//
2041	iconcount++;
2042	}
2043	}
2044	else{
2045	break;
2046	}
2047	}
2048	if (!btrigger) break;
2049	}
2050	if (btrigger){
2051
2052	//
2053	// fill the MMcEvt with all information
2054	//
2055
2056	if(!flagstoring)
2057	nstoredevents++;
2058	flagstoring = 1;
2059
2060	if (Write_McEvt) {
2061	Float_t ftime, ltime;
2062	if (reflector_file_version<6){
2063	mcevth[0].get_times(&ftime, &ltime);
2064	McEvt[0]->Fill( 0,
2065	(UShort_t) mcevth[0].get_primary() ,
2066	mcevth[0].get_energy(),
2067	-1.0,
2068	-1.0,
2069	-1.0,
2070	mcevth[0].get_theta(),
2071	mcevth[0].get_phi(),
2072	mcevth[0].get_core(),
2073	mcevth[0].get_coreX(),
2074	mcevth[0].get_coreY(),
2075	impactD[0],
2076	phiCT[0],
2077	thetaCT[0],
2078	ftime,
2079	ltime,
2080	0,
2081	0,
2082	0,
2083	0,
2084	0,
2085	0,
2086	0,
2087	(UInt_t)mcevth[0].get_CORSIKA(),
2088	(UInt_t)mcevth[0].get_AtmAbs(),
2089	(UInt_t)(mcevth[0].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()),
2090	(UInt_t) ncph[0],
2091	(UInt_t) inumphe_CT[0],
2092	(UInt_t) inumphensb[0]+inumphe_CT[0],
2093	-1.0,
2094	-1.0,
2095	-1.0);
2096	}
2097	else{
2098	Float_t Nmax, t0, tmax, a, b, c, chi2;
2099	mcevth_2[0].get_times(&ftime, &ltime);
2100	chi2=mcevth_2[0].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2101	McEvt[0]->Fill((UInt_t) mcevth_2[0].get_evt_number(),
2102	(UShort_t) mcevth_2[0].get_primary() ,
2103	mcevth_2[0].get_energy(),
2104	mcevth_2[0].get_thick0(),
2105	mcevth_2[0].get_first_target(),
2106	mcevth_2[0].get_z_first_int(),
2107	mcevth_2[0].get_theta(),
2108	mcevth_2[0].get_phi(),
2109	mcevth_2[0].get_core(),
2110	mcevth_2[0].get_coreX(),
2111	mcevth_2[0].get_coreY(),
2112	impactD[0],
2113	mcevth_2[0].get_phi_CT(),
2114	mcevth_2[0].get_theta_CT(),
2115	ftime,
2116	ltime,
2117	Nmax,
2118	t0,
2119	tmax,
2120	a,
2121	b,
2122	c,
2123	chi2,
2124	(UInt_t)mcevth_2[0].get_CORSIKA(),
2125	(UInt_t)mcevth_2[0].get_AtmAbs(),
2126	(UInt_t)(mcevth_2[0].get_MirrAbs()+mcevth_2[0].get_OutOfMirr()+mcevth_2[0].get_BlackSpot()),
2127	(UInt_t) ncph[0],
2128	(UInt_t) inumphe_CT[0],
2129	(UInt_t) inumphensb[0]+inumphe_CT[0],
2130	mcevth_2[0].get_ElecFraction(),
2131	mcevth_2[0].get_MuonFraction(),
2132	mcevth_2[0].get_OtherFraction());
2133	}
2134	}
2135	// Fill the Tree with the current leaves of each branch
2136	i=EvtTree.Fill() ;
2137
2138	// Clear the branches
2139	if(Write_McTrig){
2140	for(int i=0;i<numberBranches;i++){
2141	McTrig[i]->Clear() ;
2142	}
2143	}
2144	if( Write_RawEvt ){
2145	for(int i=0;i<numberBranches;i++){
2146	EvtHeader[i]->Clear() ;
2147	EvtData[i]->DeletePixels();
2148	}
2149	}
2150	if (Write_McEvt)
2151	McEvt[0]->Clear() ;
2152	}
2153	}
2154	// We study a single trigger condition
2155	else {
2156
2157	// Set to zero the flag to know if some conditon has triggered
2158	btrigger=0;
2159	flagstoring = 0;
2160
2161	for(int ict=0;ict<ct_Number;ict++){
2162
2163	// Setting trigger conditions
2164	Trigger_CT[ict]->SetMultiplicity(Trigger_multiplicity[ict]);
2165	Trigger_CT[ict]->SetTopology(Trigger_topology[ict]);
2166	for (int i=0;i<ct_NPixels;i++)
2167	fpixelthres[i]=qThreshold[ict][i];
2168
2169	// Rise the discrimnator threshold to avoid huge rates
2170	if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
2171	for(int ii=0;ii<ct_NPixels;ii++){
2172	if( nsb_phepns_rotated[ict][ii]>riseDiskThres)
2173	fpixelthres[ii]=secureDiskThres;
2174	}
2175
2176	Trigger_CT[ict]->SetThreshold(fpixelthres);
2177
2178	Trigger_CT[ict]->Diskriminate() ;
2179
2180	//
2181	// look if in all the signals in the trigger signal branch
2182	// is a possible Trigger. Therefore we habe to diskriminate all
2183	// the simulated analog signals (Method Diskriminate in class
2184	// MTrigger). We look simultanously for the moments at which
2185	// there are more than TRIGGER_MULTI pixels above the
2186	// CHANNEL_THRESHOLD.
2187	//
2188
2189	Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ;
2190
2191	Lev1MT[ict] = 0 ;
2192
2193	//
2194	// Start the First Level Trigger simulation
2195	//
2196
2197	if ( Lev0MT[ict] > 0 \|\| Write_All_Images) {
2198
2199	Lev1MT[ict]= Trigger_CT[ict]->FirstLevel();
2200	}
2201	if (Lev1MT[ict]>0){
2202	++ntrigger[ict];
2203	}
2204	}
2205
2206	Int_t NumImages = 0;
2207	Int_t CT_triggered=0;
2208	for(int ict=0;ict<ct_Number;ict++){
2209	if(NumImages==0 && Lev1MT[ict]>0)
2210	CT_triggered=ict;
2211	NumImages = (NumImages>=Lev1MT[ict])?NumImages:1;
2212	Lev0MT[ict]=1;
2213	if (Lev1MT[ict]==0 && Write_All_Images){
2214	NumImages=1;
2215	Lev0MT[ict]=0;
2216	}
2217	}
2218
2219	for(int ict=0;ict<ct_Number;ict++){
2220	for(Int_t ii=0;ii<NumImages;ii++){
2221
2222	btrigger=1;
2223
2224	// Loop over different level one triggers
2225
2226	//
2227	// fill inside class fadc the member output
2228	//
2229	if(Lev1MT[ict]>0)
2230	Fadc_CT[ict]->
2231	TriggeredFadc(Trigger_CT[ict]->GetFirstLevelTime(ii));
2232	else
2233	Fadc_CT[ict]->
2234	TriggeredFadc(Trigger_CT[CT_triggered]->GetFirstLevelTime(ii));
2235
2236	if(!flagstoring)
2237	nstoredevents++;
2238	flagstoring = 1;
2239
2240	if (Write_McTrig){
2241	McTrig[ict]->SetFirstLevel (Lev1MT[ict]);
2242	McTrig[ict]
2243	->SetTime(Trigger_CT[ict]->GetFirstLevelTime(ii),ii+1);
2244	Trigger_CT[ict]->GetMapDiskriminator(trigger_map);
2245	McTrig[ict]->SetMapPixels(trigger_map,ii);
2246	}
2247
2248	// Fill Evt information
2249
2250	if (Write_RawEvt){
2251
2252	//
2253	// Fill the header of this event
2254	//
2255
2256	EvtHeader[ict]
2257	->FillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ;
2258
2259	// fill pixel information
2260
2261	if (Lev1MT[ict] \|\| Write_All_Images){
2262	if (addElecNoise) Fadc_CT[ict]->DigitalNoise();
2263	for(UInt_t i=0;
2264	i<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
2265	i++){
2266	if(!Fadc_CT[ict]->IsPixelUsed(i)) continue;
2267	for (j=0;j<FADC_SLICES;j++){
2268	fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(i,j),j);
2269	fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(i,j),j);
2270	}
2271	EvtData[ict]->AddPixel(i,fadcValues,0);
2272	EvtData[ict]->AddPixel(i,fadcValuesLow,kTRUE);
2273	}
2274	}
2275	}
2276	}
2277	//
2278	// if a first level trigger occurred, then
2279	// 1. do some other stuff (not implemented)
2280	// 2. start the gui tool
2281
2282	if(FADC_Scan){
2283	if ( Lev0MT[ict] > 0 ) {
2284	Fadc_CT[ict]->ShowSignal( McEvt[ict], (Float_t) 60. ) ;
2285	}
2286	}
2287
2288	if(Trigger_Scan){
2289	if ( Lev0MT[ict] > 0 ) {
2290	Trigger_CT[ict]->ShowSignal(McEvt[ict]) ;
2291	}
2292	}
2293
2294	} // end CT loop
2295
2296	// If there is trigger in some telecope or we store all showers
2297	if(btrigger){
2298	if (Write_McEvt){
2299	//
2300	// fill the MMcEvt with all information
2301	//
2302
2303	for (int ict=0;ict<ct_Number;ict++){
2304	Float_t ftime, ltime;
2305	if (reflector_file_version<6){
2306	mcevth[ict].get_times(&ftime, &ltime);
2307	McEvt[ict]->Fill( 0,
2308	(UShort_t) mcevth[0].get_primary() ,
2309	mcevth[ict].get_energy(),
2310	-1.0,
2311	-1.0,
2312	-1.0,
2313	mcevth[ict].get_theta(),
2314	mcevth[ict].get_phi(),
2315	mcevth[ict].get_core(),
2316	mcevth[ict].get_coreX(),
2317	mcevth[ict].get_coreY(),
2318	impactD[ict],
2319	phiCT[ict],
2320	thetaCT[ict],
2321	ftime,
2322	ltime,
2323	0,
2324	0,
2325	0,
2326	0,
2327	0,
2328	0,
2329	0,
2330	(UInt_t)mcevth[ict].get_CORSIKA(),
2331	(UInt_t)mcevth[ict].get_AtmAbs(),
2332	(UInt_t)(mcevth[ict].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()),
2333	(UInt_t) ncph[ict],
2334	(UInt_t) inumphe_CT[ict],
2335	(UInt_t) inumphensb[ict]+inumphe_CT[ict],
2336	-1.0,
2337	-1.0,
2338	-1.0);
2339	}
2340	else{
2341	Float_t Nmax, t0, tmax, a, b, c, chi2;
2342	mcevth_2[ict].get_times(&ftime, &ltime);
2343	chi2=mcevth_2[ict].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2344	McEvt[ict]->Fill( (UInt_t) mcevth_2[ict].get_evt_number(),
2345	(UShort_t) mcevth_2[ict].get_primary() ,
2346	mcevth_2[ict].get_energy(),
2347	mcevth_2[ict].get_thick0(),
2348	mcevth_2[ict].get_first_target(),
2349	mcevth_2[ict].get_z_first_int(),
2350	mcevth_2[ict].get_theta(),
2351	mcevth_2[ict].get_phi(),
2352	mcevth_2[ict].get_core(),
2353	mcevth_2[ict].get_coreX(),
2354	mcevth_2[ict].get_coreY(),
2355	impactD[ict],
2356	mcevth_2[ict].get_phi_CT(),
2357	mcevth_2[ict].get_theta_CT(),
2358	ftime,
2359	ltime,
2360	Nmax,
2361	t0,
2362	tmax,
2363	a,
2364	b,
2365	c,
2366	chi2,
2367	(UInt_t)mcevth_2[ict].get_CORSIKA(),
2368	(UInt_t)mcevth_2[ict].get_AtmAbs(),
2369	(UInt_t) (mcevth_2[ict].get_MirrAbs()+mcevth_2[ict].get_OutOfMirr()+mcevth_2[ict].get_BlackSpot()),
2370	(UInt_t) ncph[ict],
2371	(UInt_t) inumphe_CT[ict],
2372	(UInt_t) inumphensb[ict]+inumphe_CT[ict],
2373	mcevth_2[ict].get_ElecFraction(),
2374	mcevth_2[ict].get_MuonFraction(),
2375	mcevth_2[ict].get_OtherFraction());
2376	}
2377	}
2378	}
2379	// We don not count photons out of the camera.
2380
2381
2382	//
2383	// write it out to the file outfile
2384	//
2385
2386	EvtTree.Fill() ;
2387	}
2388
2389	// clear all
2390	for(int ict=0;ict<ct_Number;ict++){
2391	if (Write_RawEvt) EvtHeader[ict]->Clear() ;
2392	if (Write_RawEvt) EvtData[ict]->DeletePixels();
2393	if (Write_McTrig) McTrig[ict]->Clear() ;
2394	if (Write_McEvt) McEvt[ict]->Clear() ;
2395	}
2396	}
2397
2398	#ifdef __DEBUG__
2399	printf("\n");
2400
2401	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
2402
2403	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
2404
2405	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
2406
2407	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
2408
2409	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
2410	(int)pixels[ici][icj][PIXNUM],
2411	pixels[ici][icj][PIXX],
2412	pixels[ici][icj][PIXY],
2413	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
2414
2415	}
2416	}
2417	}
2418	}
2419
2420	for (int i=0;
2421	i<((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels(); ++i) {
2422	printf("%d (%d): ", i, npixneig[i]);
2423	for (j=0; j<npixneig[i]; ++i)
2424	printf(" %d", pixneig[i][j]);
2425	printf("\n");
2426	}
2427
2428	#endif // __DEBUG__
2429
2430
2431	// We search the maximum impact parameter fo the simualted showers
2432	maxpimpact=maxpimpact<impactD[0]?impactD[0]:maxpimpact;
2433
2434	// look for the next event
2435
2436	for(int ict=0;ict<ct_Number;ict++)
2437	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2438
2439	} // end while there is a next event
2440
2441	if( !isA( flag, FLAG_END_OF_RUN )){
2442	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
2443	break;
2444	}
2445	else { // found end of run
2446	ntshow += nshow;
2447	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
2448
2449	//fread( flag, SIZE_OF_FLAGS, 1, inputfile );
2450	for(int ict=0;ict<ct_Number;ict++)
2451	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2452
2453	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
2454	log(SIGNATURE, "End of file . . .\n");
2455	still_in_loop = FALSE;
2456	log(SIGNATURE, "Reading ASCII files at the end of the reflector file. . .\n");
2457	for(int ict=0;ict<ct_Number;ict++){
2458	read_ascii(inputfile[ict], McConfigRunHeader[ict]);
2459	}
2460	if ((! Data_From_STDIN) && ( !feof(inputfile[0]) )){
2461
2462	// we have concatenated input files.
2463	// get signature of the next part and check it.
2464
2465	if((reflector_file_version=check_reflector_file( inputfile[0] ))==FALSE){
2466	log(SIGNATURE, "Next file is not recognised as a reflector file.\n");
2467	log(SIGNATURE, "Stopping ...\n");
2468	break;
2469	}
2470
2471	}
2472
2473	for(int ict=0;ict<ct_Number;ict++)
2474	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2475
2476	} // end if found end of file
2477	} // end if found end of run
2478
2479	} // end if else found start of run
2480	} // end big while loop
2481
2482	//<@ Finally we should fill th McRunHeader
2483
2484	Float_t heights[10];
2485	time_t ltime;
2486	Float_t ftime;
2487	Float_t rnum;
2488	Float_t viewcone[2]={0,0};
2489
2490	get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
2491	if (reflector_file_version<6){
2492	telestheta=-10.0;
2493	telesphi=-10.0;
2494	mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2495	mcevth[0].get_phi_range(&shphimin,&shphimax);
2496	mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2497	mcevth[0].get_phi_range(&shphimin,&shphimax);
2498	corsika=UInt_t(mcevth[0].get_VersionPGM()*1000);
2499	for (int i=0; i< 10;i++)
2500	heights[i]=mcevth[0].get_HeightLev (i);
2501	rnum=mcevth[0].get_RunNumber();
2502	}
2503	else{
2504	telestheta=mcevth_2[0].get_theta_CT();
2505	telesphi=mcevth_2[0].get_phi_CT();
2506	mcevth_2[0].get_theta_range(&shthetamin, &shthetamax);
2507	mcevth_2[0].get_phi_range(&shphimin,&shphimax);
2508	corsika=UInt_t(mcevth_2[0].get_VersionPGM()*1000);
2509	for (int i=0; i< 10;i++)
2510	heights[i]=mcevth_2[0].get_HeightLev (i);
2511	rnum=mcevth_2[0].get_RunNumber();
2512	mcevth_2[0].get_viewcone(&viewcone[0],&viewcone[1]);
2513	}
2514
2515	get_source_off(&sofftheta,&soffphi);
2516	if(!Trigger_Loop) icontrigger=0;
2517	time (&ltime);
2518	ftime = ((Float_t)ltime)/1000;
2519
2520	if (reflector_file_version<6)
2521	McRunHeader->Fill(rnum,
2522	(UInt_t) 0,
2523	mcevth[0].get_DateRun(),
2524	ftime,
2525	icontrigger,
2526	!Write_All_Images,
2527	Write_McEvt,
2528	Write_McTrig,
2529	Write_McFADC,
2530	Write_RawEvt,
2531	addElecNoise,
2532	ct_NPixels,
2533	(UInt_t)ntshow,
2534	(UInt_t)nstoredevents,
2535	0,
2536	sfRaH,
2537	sfRaM,
2538	sfRaS,
2539	sfDeD,
2540	sfDeM,
2541	sfDeS,
2542	meanNSB,
2543	telestheta,
2544	telesphi,
2545	sofftheta,
2546	soffphi,
2547	shthetamax,
2548	shthetamin,
2549	shphimax,
2550	shphimin,
2551	maxpimpact,
2552	mcevth[0].get_CWaveLower(),
2553	mcevth[0].get_CWaveUpper(),
2554	mcevth[0].get_slope(),
2555	1,
2556	heights,
2557	corsika,
2558	(UInt_t)(reflector_file_version*100),
2559	(UInt_t)(VERSION*100),
2560	0);
2561	else
2562	McRunHeader->Fill(rnum,
2563	(UInt_t) 0,
2564	mcevth_2[0].get_DateRun(),
2565	ftime,
2566	icontrigger,
2567	!Write_All_Images,
2568	Write_McEvt,
2569	Write_McTrig,
2570	Write_McFADC,
2571	Write_RawEvt,
2572	addElecNoise,
2573	ct_NPixels,
2574	(UInt_t)ntshow,
2575	(UInt_t)nstoredevents,
2576	0,
2577	sfRaH,
2578	sfRaM,
2579	sfRaS,
2580	sfDeD,
2581	sfDeM,
2582	sfDeS,
2583	meanNSB,
2584	telestheta,
2585	telesphi,
2586	sofftheta,
2587	soffphi,
2588	shthetamax,
2589	shthetamin,
2590	shphimax,
2591	shphimin,
2592	maxpimpact,
2593	mcevth_2[0].get_CWaveLower(),
2594	mcevth_2[0].get_CWaveUpper(),
2595	mcevth_2[0].get_slope(),
2596	1,
2597	heights,
2598	corsika,
2599	(UInt_t)(reflector_file_version*100),
2600	(UInt_t)(VERSION*100),
2601	0);
2602	// Fill some missing values for MRawRunHeader
2603
2604	RunHeader->SetRunNumber(rnum);
2605	RunHeader->SetNumEvents(nstoredevents);
2606
2607	// Fill MMcCorsikaRunHeader
2608
2609	Float_t constantC[50];
2610	mcrunh.get_constantC(&constantC[0]);
2611	Float_t constantCKA[40];
2612	mcrunh.get_constantCKA(&constantCKA[0]);
2613	Float_t constantCETA[5];
2614	mcrunh.get_constantCETA(&constantCETA[0]);
2615	Float_t constantCSTRBA[11];
2616	mcrunh.get_constantCSTRBA(&constantCSTRBA[0]);
2617	Float_t constantAATM[5];
2618	mcrunh.get_constantAATM(&constantAATM[0]);
2619	Float_t constantBATM[5];
2620	mcrunh.get_constantBATM(&constantBATM[0]);
2621	Float_t constantCATM[5];
2622	mcrunh.get_constantCATM(&constantCATM[0]);
2623	Float_t constantNFL[4];
2624	mcrunh.get_constantNFL(&constantNFL[0]);
2625
2626	if(reflector_file_version>5)
2627	McCorsikaRunHeader->Fill(rnum,
2628	mcrunh.get_date(),
2629	corsika,
2630	1,
2631	heights,
2632	mcevth_2[0].get_slope(),
2633	mcrunh.get_ELow(),
2634	mcrunh.get_EUpp(),
2635	mcrunh.get_EGS4(),
2636	mcrunh.get_NKG(),
2637	mcrunh.get_Ecutoffh(),
2638	mcrunh.get_Ecutoffm(),
2639	mcrunh.get_Ecutoffe(),
2640	mcrunh.get_Ecutoffg(),
2641	constantC,
2642	constantCKA,
2643	constantCETA,
2644	constantCSTRBA,
2645	constantAATM,
2646	constantBATM,
2647	constantCATM,
2648	constantNFL,
2649	viewcone,
2650	mcrunh.get_wobble(),
2651	mcrunh.get_atmophere()
2652	);
2653
2654	// Store qe for each PMT in output file
2655	TArrayF qe_pmt;
2656	TArrayF wav_pmt;
2657
2658	for(int ict=0;ict<ct_Number;ict++){
2659	McConfigRunHeader[ict]->InitSizePMTs(ct_NPixels);
2660	for(int i=0; i<(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();i++){
2661	McConfigRunHeader[ict]->AddPMT(i);
2662	MGeomPMT &pmt = McConfigRunHeader[ict]->GetPMT(i);
2663	qe_pmt.Set(pointsQE[ict],QE[ict][i][1]);
2664	wav_pmt.Set(pointsQE[ict],QElambda);
2665	pmt.SetArraySize(pointsQE[ict]);
2666	pmt.SetPMTContent(i,wav_pmt,qe_pmt);
2667	}
2668
2669	// Store Light Guides factors in the output file
2670	TArrayF theta_lg;
2671	TArrayF factor_lg;
2672	theta_lg.Set(pointsWC,WC[0]);
2673	factor_lg.Set(pointsWC,WC[1]);
2674	McConfigRunHeader[ict]->SetLightGuides(theta_lg,factor_lg);
2675
2676	}
2677
2678	// Fill the Header Tree with the current leaves of each branch
2679	HeaderTree.Fill() ;
2680
2681	//++
2682	// put the Event to the root file
2683	//--
2684
2685	outfile_temp.Write() ;
2686	outfile_temp.Close() ;
2687
2688	// close input file
2689
2690	if (Trigger_Loop){
2691	log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
2692	ntshow, ntcph[0] );
2693	datafile<<ntshow<<" event(s), with a total of "<<ntcph[0]<<" C.photons"<<endl;
2694	log( SIGNATURE, "Trigger Mode. \n");
2695	log( SIGNATURE, "Fraction of triggers: \n");
2696	datafile<<"Fraction of triggers: "<<endl;
2697	for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
2698	for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
2699	for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
2700	log( SIGNATURE, "Thres %5.1f, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
2701	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);
2702	datafile<<"Thres "<<fthrescount<<", Multi "<<imulticount<<", Topo"<<isorttopo[itopocount]<<": ";
2703	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;
2704	}
2705	}
2706	}
2707	}
2708	else{
2709	for(int ict=0;ict<ct_Number;ict++){
2710	log( SIGNATURE,
2711	"%d event(s), with a total of %d C.photons in CT %i (%s)\n",
2712	ntshow, ntcph[ict],ict,GeometryName[ict] );
2713	datafile<<ntshow<<" event(s), with a total of "<<ntcph[ict]
2714	<<" C.photons in CT "<<ict<<" ("<<GeometryName[ict]<<")"<<endl;
2715	log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
2716	((float)ntrigger[ict]) / ((float)ntshow) * 100.0, ntrigger[ict], ntshow);
2717	datafile<<"Fraction of triggers: "<<((float)ntrigger[ict]) / ((float)ntshow) * 100.0<<" ("<<ntrigger[ict]<<" out of "<<ntshow<<" )"<<endl;
2718	}
2719	}
2720
2721	// close files
2722
2723	log( SIGNATURE, "Closing files\n" );
2724
2725	if( ! Data_From_STDIN ){
2726	for(int ict=0;ict<ct_Number;ict++)
2727	fclose( inputfile[ict] );
2728	}
2729	datafile.close();
2730
2731	// program finished
2732
2733	log( SIGNATURE, "Done.\n");
2734
2735	return( 0 );
2736	}
2737	//!@}
2738
2739	// @T \newpage
2740
2741	//!@subsection Functions definition.
2742
2743	//!-----------------------------------------------------------
2744	// @name present
2745	//
2746	// @desc Make some presentation
2747	//
2748	// @date Sat Jun 27 05:58:56 MET DST 1998
2749	//------------------------------------------------------------
2750	// @function
2751
2752	//!@{
2753	void
2754	present(void)
2755	{
2756	cout << "##################################################\n"
2757	<< SIGNATURE << '\n' << '\n'
2758	<< "Processor of the reflector output\n"
2759	<< "J C Gonzalez, Jun 1998\n"
2760	<< "##################################################\n\n"
2761	<< flush ;
2762	}
2763	//!@}
2764
2765
2766	//!-----------------------------------------------------------
2767	// @name usage
2768	//
2769	// @desc show help
2770	//
2771	// @date Tue Dec 15 16:23:30 MET 1998
2772	//------------------------------------------------------------
2773	// @function
2774
2775	//!@{
2776	void
2777	usage(void)
2778	{
2779	present();
2780	cout << "\nusage ::\n\n"
2781	<< "\t camera "
2782	<< " [ -@ paramfile ] "
2783	<< " [ -h ] "
2784	<< "\n\n or \n\n"
2785	<< "\t camera < paramfile"
2786	<< "\n\n";
2787	exit(0);
2788	}
2789	//!@}
2790
2791
2792	//!-----------------------------------------------------------
2793	// @name log
2794	//
2795	// @desc function to send log information
2796	//
2797	// @var funct Name of the caller function
2798	// @var fmt Format to be used (message)
2799	// @var ... Other information to be shown
2800	//
2801	// @date Sat Jun 27 05:58:56 MET DST 1998
2802	//------------------------------------------------------------
2803	// @function
2804
2805	//!@{
2806	void
2807	log(const char funct, char fmt, ...)
2808	{
2809	va_list args;
2810
2811	// Display the name of the function that called error
2812	printf("[%s]: ", funct);
2813
2814	// Display the remainder of the message
2815	va_start(args, fmt);
2816	vprintf(fmt, args);
2817	va_end(args);
2818	}
2819	//!@}
2820
2821
2822	//!-----------------------------------------------------------
2823	// @name error
2824	//
2825	// @desc function to send an error message, and abort the program
2826	//
2827	// @var funct Name of the caller function
2828	// @var fmt Format to be used (message)
2829	// @var ... Other information to be shown
2830	//
2831	// @date Sat Jun 27 05:58:56 MET DST 1998
2832	//------------------------------------------------------------
2833	// @function
2834
2835	//!@{
2836	void
2837	error(const char funct, char fmt, ...)
2838	{
2839	va_list args;
2840
2841	// Display the name of the function that called error
2842	fprintf(stdout, "ERROR in %s: ", funct);
2843
2844	// Display the remainder of the message
2845	va_start(args, fmt);
2846	vfprintf(stdout, fmt, args);
2847	va_end(args);
2848
2849	perror(funct);
2850
2851	exit(1);
2852	}
2853	//!@}
2854
2855
2856	//!-----------------------------------------------------------
2857	// @name isA
2858	//
2859	// @desc returns TRUE(FALSE), if the flag is(is not) the given
2860	//
2861	// @var s1 String to be searched
2862	// @var flag Flag to compare with string s1
2863	// @return TRUE: both strings match; FALSE: oth.
2864	//
2865	// @date Wed Jul 8 15:25:39 MET DST 1998
2866	//------------------------------------------------------------
2867	// @function
2868
2869	//!@{
2870	int
2871	isA( char * s1, const char * flag ) {
2872	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
2873	}
2874	//!@}
2875
2876
2877	//!-----------------------------------------------------------
2878	// @name read_QE
2879	//
2880	// @desc read QE data
2881	//
2882	// @date thu 5 17:59:57 CEST 2002
2883	//------------------------------------------------------------
2884	// @function
2885
2886	//!@{
2887	void
2888	read_QE(char fname[256], int ict){
2889	ifstream qefile;
2890	char line[LINE_MAX_LENGTH];
2891	int i, j, icount;
2892	float qe;
2893
2894	//------------------------------------------------------------
2895	// second, pixels' QE
2896
2897	// try to open the file
2898
2899	log("read_QE", "Opening the file \"%s\" . . .\n", fname);
2900
2901	qefile.open( fname );
2902
2903	// if it is wrong or does not exist, exit
2904
2905	if ( qefile.bad() )
2906	error( "read_QE", "Cannot open \"%s\". Exiting.\n", fname );
2907
2908	// read file
2909
2910	log("read_QE", "Reading data . . .\n");
2911
2912	i=-1;
2913	icount = 0;
2914
2915	while ( ! qefile.eof() ) {
2916
2917	// get line from the file
2918
2919	qefile.getline(line, LINE_MAX_LENGTH);
2920
2921	// skip if comment
2922
2923	if ( *line == '#' )
2924	continue;
2925
2926	// if it is the first valid value, it is the number of QE data points
2927
2928	if ( i < 0 ) {
2929
2930	// get the number of datapoints
2931
2932	sscanf(line, "%d", &pointsQE[ict]);
2933
2934	// allocate memory for the table of QEs
2935
2936	QE[ict] = new float ** [ct_NPixels];
2937
2938	for ( i=0; i<ct_NPixels; ++i ) {
2939	QE[ict][i] = new float * [2];
2940	QE[ict][i][0] = new float[pointsQE[ict]];
2941	QE[ict][i][1] = new float[pointsQE[ict]];
2942	}
2943
2944	QElambda = new float [pointsQE[ict]];
2945
2946	for ( i=0; i<pointsQE[ict]; ++i ) {
2947	qefile.getline(line, LINE_MAX_LENGTH);
2948	sscanf(line, "%f", &QElambda[i]);
2949	}
2950
2951	i=0;
2952
2953	continue;
2954	}
2955
2956	// get the values (num-pixel, num-datapoint, QE-value)
2957
2958	if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2959	break;
2960
2961	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2962	((j-1) < pointsQE[ict]) && ((j-1) > -1) ) {
2963	QE[ict][i-1][0][j-1] = QElambda[j-1];
2964	QE[ict][i-1][1][j-1] = qe;
2965	}
2966
2967	if ( i > ct_NPixels) break;
2968
2969	icount++;
2970
2971	}
2972
2973	if(icount/pointsQE[ict] < ct_NPixels){
2974	error( "read_QE", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2975	icount/pointsQE[ict], ct_NPixels );
2976	}
2977
2978	// close file
2979
2980	qefile.close();
2981
2982	// test QE
2983
2984	for(icount=0; icount< ct_NPixels; icount++){
2985	for(i=0; i<pointsQE[ict]; i++){
2986	if( QE[ict][icount][0][i] < 100. \|\| QE[ict][icount][0][i] > 1000. \|\|
2987	QE[ict][icount][1][i] < 0. \|\| QE[ict][icount][1][i] > 100.){
2988	error( "read_QE", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2989	icount, i, QE[ict][icount][0][i], QE[ict][icount][1][i] );
2990	}
2991	}
2992	}
2993
2994	// end
2995
2996	log("read_QE", "Done.\n");
2997	}
2998	//!@}
2999
3000	//!-----------------------------------------------------------
3001	// @name read_WC
3002	//
3003	// @desc read WC data
3004	//
3005	// @date thu 5 17:59:57 CEST 2002
3006	//------------------------------------------------------------
3007	// @function
3008
3009	//!@{
3010	void
3011	read_WC(void){
3012	ifstream wcfile;
3013	char line[LINE_MAX_LENGTH];
3014	int i;
3015
3016	//------------------------------------------------------------
3017	// Read Light Guides data
3018
3019	// try to open the file
3020
3021	log("read_WC", "Opening the file \"%s\" . . .\n", WC_FILE);
3022
3023	wcfile.open( WC_FILE );
3024
3025	// if it is wrong or does not exist, exit
3026
3027	if ( wcfile.bad() )
3028	error( "read_WC", "Cannot open \"%s\". Exiting.\n", WC_FILE );
3029
3030	// read file
3031
3032	log("read_WC", "Reading data . . .\n");
3033
3034	// get line from the file
3035
3036	wcfile.getline(line, LINE_MAX_LENGTH);
3037
3038	// get the number of datapoints
3039
3040	sscanf(line, "%d", &pointsWC);
3041
3042	// allocate memory for the table of QEs
3043
3044	WC = new float * [2];
3045	WC[0] = new float[pointsWC];
3046	WC[1] = new float[pointsWC];
3047
3048	for ( i=0; i<pointsWC; ++i ) {
3049	wcfile.getline(line, LINE_MAX_LENGTH);
3050	sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
3051	}
3052
3053	// close file
3054
3055	wcfile.close();
3056
3057	// read
3058
3059	log("read_WC", "Done.\n");
3060	}
3061	//!@}
3062
3063
3064	//!-----------------------------------------------------------
3065	// @name read_ascii
3066	//
3067	// @desc read ascii configuration files used by the reflector
3068	//
3069	// @date tue dec 10 17:14:10 CET 2002
3070	//------------------------------------------------------------
3071	// @function
3072
3073	//!@{
3074	void
3075	read_ascii(FILE sp, MMcConfigRunHeader config)
3076	{
3077	Float_t radius = -1.0;
3078	Float_t focal = -1.0;
3079	Float_t stdfocal = -1.0;
3080	Float_t point = -1.0;
3081	Float_t stdpoint = -1.0;
3082	Float_t adjust = -1.0;
3083	Float_t spot = -1.0;
3084	Float_t camwidth = -1.0;
3085
3086	Int_t imir;
3087	Float_t f,sx,sy,x,y,z,thetan,phin,xn,yn,zn;
3088	Float_t dx,dy;
3089	Int_t nummir, numref;
3090	Float_t wav,ref;
3091
3092	Char_t token[40];
3093	Char_t line[511];
3094	Char_t flag;
3095
3096	while(1){
3097	if((flag=fgetc(sp))==EOF)
3098	break;
3099	fgets(&line[1],500,sp);
3100	line[0]=flag;
3101	if (line[0]== '#' && line[2]=='n' && line[3]=='u' && line[4]=='m' &&
3102	line[5]== 'b' && line[12]=='d' && line[13]=='a' && line[14]=='t'){
3103	fgets(line,500,sp);
3104	sscanf(line,"%i",&numref);
3105	TArrayF wavarray(numref);
3106	TArrayF refarray(numref);
3107	fgets(line,500,sp);
3108	for(int i=0; i<numref;i++){
3109	fgets(line,500,sp);
3110	sscanf(line,"%f %f",&wav,&ref);
3111	wavarray[i]=wav;
3112	refarray[i]=ref;
3113	}
3114	for (int j=0; j<nummir;j++){
3115	MGeomMirror &mirror = config->GetMirror(j);
3116	mirror.SetArraySize(numref);
3117	mirror.SetReflectivity(wavarray, refarray);
3118	}
3119	continue;
3120	}
3121	if (line[0]== '#')
3122	continue;
3123	if (strstr(line, "type") == line)
3124	continue;
3125	if (strstr(line, "focal_distance") == line){
3126	sscanf(line,"%s %f",token,&focal);
3127	continue;
3128	}
3129	if (strstr(line, "focal_std") == line){
3130	sscanf(line,"%s %f",token,&stdfocal);
3131	continue;
3132	}
3133	if (strstr(line, "point_spread") == line){
3134	sscanf(line,"%s %f",token,&point);
3135	continue;
3136	}
3137	if (strstr(line, "point_std") == line){
3138	sscanf(line,"%s %f",token,&stdpoint);
3139	continue;
3140	}
3141	if (strstr(line, "adjustment_dev") == line){
3142	sscanf(line,"%s %f",token,&adjust);
3143	continue;
3144	}
3145	if (strstr(line, "black_spot") == line){
3146	sscanf(line,"%s %f",token,&spot);
3147	continue;
3148	}
3149	if (strstr(line, "camera_width") == line){
3150	sscanf(line,"%s %f",token,&camwidth);
3151	continue;
3152	}
3153	if (strstr(line, "n_pixels") == line)
3154	continue;
3155	if (strstr(line, "pixel_width") == line)
3156	continue;
3157	if (strstr(line, "n_centralpixels") == line)
3158	continue;
3159	if (strstr(line, "n_gappixels") == line)
3160	continue;
3161	if (strstr(line, "n_mirrors") == line){
3162	sscanf(line,"%s %i",token,&nummir);
3163	config->InitSizeMirror(nummir);
3164	continue;
3165	}
3166	if (strstr(line, "r_mirror") == line){
3167	sscanf(line,"%s %f",token,&radius);
3168	continue;
3169	}
3170	if (strstr(line, "define_mirrors") == line){
3171	for(int i=0;i<nummir;i++){
3172	fgets(line,500,sp);
3173	sscanf(line,"%i %f %f %f %f %f %f %f %f %f %f %f",
3174	&imir,&f,&sx,&sy,&x,&y,&z,&thetan,&phin,&xn,&yn,&zn);
3175	config->AddMirror(i);
3176	MGeomMirror &mirror = config->GetMirror(i);
3177	mirror.SetMirrorContent(imir,f,sx,sy,x,y,z,thetan,phin,xn,yn,zn);
3178	}
3179	fgets(line,500,sp);
3180	while(line[2]!='w' \|\| line[3]!='i' \|\| line[4]!='t' \|\| line[5]!='h'){
3181	fgets(line,500,sp);
3182	}
3183	fgets(line,500,sp);
3184	for(int i=0;i<nummir;i++){
3185	fgets(line,500,sp);
3186	sscanf(line,"%f %f",&dx,&dy);
3187	MGeomMirror &mirror = config->GetMirror(i);
3188	mirror.SetMirrorDeviations(dx,dy);
3189	}
3190	continue;
3191	}
3192	}
3193	config->SetMagicDef(radius, focal, stdfocal, point,
3194	stdpoint, adjust, spot, camwidth);
3195	}
3196
3197
3198	//!-----------------------------------------------------------
3199	// @name igen_pixel_coordinates
3200	//
3201	// @desc generate the pixel center coordinates
3202	//
3203	// @var *pcam structure camera containing all the
3204	// camera information
3205	// @return total number of pixels
3206	//
3207	// DP
3208	//
3209	// @date Thu Oct 14 10:41:03 CEST 1999
3210	//------------------------------------------------------------
3211	// @function
3212
3213	//!@{
3214	/****** igen_pixel_coordinates() *******************************/
3215
3216	int igen_pixel_coordinates(struct camera *pcam) {
3217	/* generate pixel coordinates, return value is number of pixels */
3218
3219	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
3220	float fsegment_fract;
3221	double dtsize;
3222	double dhsize;
3223	double dpsize;
3224	double dxfirst_pix;
3225	double dyfirst_pix;
3226	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5;
3227	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5;
3228
3229
3230	double dstartx, dstarty; /* for the gap pixels and outer pixels */
3231	int j, nrow;
3232
3233	dpsize = pcam->dpixdiameter_cm;
3234	dtsize = dpsize * sqrt(3.) / 2.;
3235	dhsize = dpsize / 2.;
3236
3237	/* Loop over central pixels to generate co-ordinates */
3238
3239	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
3240
3241	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
3242
3243	pcam->dpixsizefactor[ipixno-1] = 1.;
3244
3245	in = 0;
3246
3247	i = 0;
3248	itot_inside_ring = 0;
3249	iring_no = 0;
3250
3251	/* Calculate the number of pixels out to and including the ring containing pixel number */
3252	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
3253
3254	while (itot_inside_ring == 0){
3255
3256	iN = 3(i(i+1)) + 1;
3257
3258	if (ipixno <= iN){
3259	iring_no = i;
3260	itot_inside_ring = iN;
3261	}
3262
3263	i++;
3264	}
3265
3266
3267	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
3268
3269	ipix_in_ring = 0;
3270	for (i = 0; i < iring_no; ++i){
3271
3272	ipix_in_ring = ipix_in_ring + 6;
3273	}
3274
3275	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
3276	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
3277	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
3278
3279	isegment = 0;
3280	fsegment_fract = 0.;
3281	if (iring_no > 0) {
3282
3283	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
3284
3285	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
3286
3287	}
3288
3289	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
3290	/* pixel width (flat to flat) dpsize. */
3291
3292	dxfirst_pix = dpsize*iring_no;
3293	dyfirst_pix = 0.;
3294
3295	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
3296
3297	ddxseg1 = - dhsize*iring_no;
3298	ddyseg1 = dtsize*iring_no;
3299	ddxseg2 = -dpsize*iring_no;
3300	ddyseg2 = 0.;
3301	ddxseg3 = ddxseg1;
3302	ddyseg3 = -ddyseg1;
3303	ddxseg4 = -ddxseg1;
3304	ddyseg4 = -ddyseg1;
3305	ddxseg5 = -ddxseg2;
3306	ddyseg5 = 0.;
3307
3308	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
3309	/* anti-clockwise around the ring by adding the segment to segment vectors. */
3310
3311	switch (isegment) {
3312
3313	case 0:
3314
3315	pcam->dxc[ipixno-1] = 0.;
3316	pcam->dyc[ipixno-1] = 0.;
3317
3318	case 1:
3319	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
3320	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
3321
3322	break;
3323
3324	case 2:
3325
3326	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
3327	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
3328
3329	break;
3330
3331	case 3:
3332
3333	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
3334	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
3335
3336	break;
3337
3338	case 4:
3339
3340	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
3341	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
3342
3343	break;
3344
3345	case 5:
3346
3347	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
3348	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
3349
3350	break;
3351
3352	case 6:
3353
3354	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
3355	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
3356
3357	break;
3358
3359	default:
3360
3361	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
3362	return(0);
3363
3364	} /* end switch */
3365
3366	} /* end for */
3367
3368	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
3369	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
3370
3371	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
3372
3373	j = pcam->inumcentralpixels;
3374
3375	for(i=0; i<pcam->inumgappixels; i=i+6){
3376	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
3377	pcam->dyc[j + i ] = dstarty;
3378	pcam->dpixsizefactor[j + i] = 1.;
3379	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
3380	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
3381	pcam->dpixsizefactor[j + i + 1] = 1.;
3382	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
3383	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
3384	pcam->dpixsizefactor[j + i+ 2] = 1.;
3385	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
3386	pcam->dyc[j + i + 3] = dstarty;
3387	pcam->dpixsizefactor[j + i+ 3] = 1.;
3388	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
3389	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
3390	pcam->dpixsizefactor[j + i+ 4] = 1.;
3391	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
3392	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
3393	pcam->dpixsizefactor[j + i + 5] = 1.;
3394	} /* end for */
3395	} /* end if */
3396
3397	/* generate positions of the outer pixels */
3398
3399	if( pcam->inumbigpixels > 0 ){
3400
3401	j = pcam->inumcentralpixels + pcam->inumgappixels;
3402
3403	for(i=0; i<pcam->inumbigpixels; i++){
3404	pcam->dpixsizefactor[j + i] = 2.;
3405	}
3406
3407	in = 0;
3408
3409	nrow = (int) ceil(dstartx / 2. / dpsize);
3410
3411	while(in < pcam->inumbigpixels){
3412
3413	pcam->dxc[j + in] = dstartx + dpsize;
3414	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
3415	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
3416	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
3417	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
3418	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
3419	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
3420	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
3421	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
3422	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
3423	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
3424	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
3425	for(i=1; i<nrow; i++){
3426	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
3427	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
3428	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
3429	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
3430	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
3431	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
3432	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
3433	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
3434	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
3435	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
3436	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
3437	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
3438	}
3439	in = in + 6 * nrow;
3440	dstartx = dstartx + 2. * dpsize;
3441	nrow = nrow + 1;
3442	} /* end while */
3443
3444	} /* end if */
3445
3446	return(pcam->inumpixels);
3447
3448	}
3449	//!@}
3450
3451	//!-----------------------------------------------------------
3452	// @name bpoint_is_in_pix
3453	//
3454	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
3455	//
3456	// @var *pcam structure camera containing all the
3457	// camera information
3458	// @var dx, dy point coordinates in centimeters
3459	// @var ipixnum pixel number (starting at 0)
3460	// @return TRUE if the point is inside the pixel, FALSE otherwise
3461	//
3462	// DP
3463	//
3464	// @date Thu Oct 14 16:59:04 CEST 1999
3465	//------------------------------------------------------------
3466	// @function
3467
3468	//!@{
3469
3470	/****** bpoint_is_in_pix() *************************************/
3471
3472	#define sqrt13 0.577350269 // = 1./sqrt(3.)
3473	#define sqrt3 1.732050807 // = sqrt(3.)
3474
3475	int bpoint_is_in_pix(double dx, double dy, MGeomCam *pgeomcam)
3476	{
3477	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
3478	/* the pixel is assumed to be a "closed set" */
3479
3480	/*
3481	a = length of one of the edges of one pixel,
3482	b = half the width of one pixel
3483	*/
3484
3485	const int numN = pgeomcam->GetNumPixels();
3486
3487	for (int i=0; i<numN; i++)
3488	{
3489	MGeomPix &pixel = (*pgeomcam)[i];
3490	const double b = pixel.GetD()/2;
3491	const double a = pixel.GetD()/sqrt3;
3492
3493	const double xx = dx - pixel.GetX();
3494	const double yy = dy - pixel.GetY();
3495
3496	if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) \|\|
3497	((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){
3498
3499	return i; // inside i
3500	}
3501
3502	// return -1; // outside
3503	}
3504
3505	return -1; // outside
3506	}
3507
3508	//!@}
3509
3510	//------------------------------------------------------------
3511	// @name dist_r_P
3512	//
3513	// @desc distance straight line r - point P
3514	//
3515	// @date Sat Jun 27 05:58:56 MET DST 1998
3516	// @function @code
3517	//------------------------------------------------------------
3518	// dist_r_P
3519	//
3520	// distance straight line r - point P
3521	//------------------------------------------------------------
3522
3523	float
3524	dist_r_P(float a, float b, float c,
3525	float l, float m, float n,
3526	float x, float y, float z)
3527	{
3528	return (
3529	sqrt((SQR((a-x)m-(b-y)l) +
3530	SQR((b-y)n-(c-z)m) +
3531	SQR((c-z)l-(a-x)n))/
3532	(SQR(l)+SQR(m)+SQR(n))
3533	)
3534	);
3535	}
3536
3537	//------------------------------------------------------------
3538	// @name check_reflector_file
3539	//
3540	// @desc check if a given reflector file has the right signature
3541	// @desc return TRUE or FALSE
3542	//
3543	// @date Mon Feb 14 16:44:21 CET 2000
3544	// @function @code
3545	//------------------------------------------------------------
3546
3547	int check_reflector_file(FILE *infile){
3548
3549	char sign[20]; // auxiliary variable
3550
3551	strcpy(sign, REFL_SIGNATURE_B);
3552
3553	fread( (char *)sign, strlen(REFL_SIGNATURE_B), 1, infile);
3554	if (strcmp(sign, REFL_SIGNATURE_A) == 0){
3555	fread( (char *)sign, 1, 1, infile);
3556	return 4;
3557	}
3558	else if (strcmp(sign, REFL_SIGNATURE_B) == 0){
3559	fread( (char *)sign, 1, 1, infile);
3560	return 5;
3561	}
3562	else if (strcmp(sign, REFL_SIGNATURE_C) == 0){
3563	// An empty bin has been removed and therefore we do not need to rezd it.
3564	return 6;
3565	}
3566	else {
3567	cout << "ERROR: Signature of .rfl file is not correct\n";
3568	cout << '"' << sign << '"' << '\n';
3569	cout << "should be: " << REFL_SIGNATURE_A <<" or "<< REFL_SIGNATURE_B <<" or " << REFL_SIGNATURE_C <<" or "<< '\n';
3570	return(FALSE);
3571	}
3572
3573
3574	}
3575
3576	//------------------------------------------------------------
3577	// @name lin_interpol
3578	//
3579	// @desc interpolate linearly between two points returning the
3580	// @desc y-value of the result
3581	//
3582	// @date Thu Feb 17 11:31:32 CET 2000
3583	// @function @code
3584	//------------------------------------------------------------
3585
3586	float lin_interpol( float x1, float y1, float x2, float y2, float x){
3587
3588	if( (x2 - x1)<=0. ){ // avoid division by zero, return average
3589	cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
3590	return((y1+y2)/2.);
3591	}
3592	else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
3593	return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
3594	}
3595	}
3596
3597	//------------------------------------------------------------
3598	// @name size_rotated
3599	//
3600	// @desc It rotates the NSB
3601	//
3602	// @date Tue Jan 7 17:09:25 CET 2003
3603	// @function @code
3604	//------------------------------------------------------------
3605
3606	int size_rotated(
3607	float *rotated,
3608	float nsb[iMAXNUMPIX],
3609	float rho)
3610	{
3611	int r=0;
3612	float size_rotated;
3613	float Num_Pixels;
3614	float PixNum=0;
3615	float rho_pixel;
3616	int j=0;
3617	int k=0;
3618
3619	for(int i=1; i<iMAXNUMPIX;i++){
3620	// Substituir per Int_t radius[iMAXNUMPIX]={0,1,1,1,1,1,1,2,2,2,2, ...}
3621	Num_Pixels=0;
3622	for (int ii=1; ii<17;ii++){
3623	if (Num_Pixels >= i){
3624	r=ii-1;
3625	break;
3626	}
3627
3628
3629
3630	if (ii<12){
3631	Num_Pixels+=ii*6;
3632	PixNum=6*ii;
3633	// size_rotated = (360/PixNum);
3634	//rho_pixel = rho/size_rotated;
3635	}
3636
3637
3638	else
3639	{
3640	Num_Pixels+=(ii-6)*6;
3641	PixNum=(ii-6)*6;
3642	// size_rotated = (360/PixNum);
3643	// rho_pixel = rho/size_rotated;
3644	}
3645	}
3646
3647
3648
3649	size_rotated = (360/PixNum);
3650	rho_pixel = rho/size_rotated;
3651
3652	//Buscar j i k que no canvin de r!!!
3653
3654	j=i-int(rho_pixel);
3655
3656	//cout<<"j inicial "<<j<<endl;
3657	int MINPIXinner[13]= {0,1,7,19,37,61,91,127,169,217,271,331,397};
3658	int MINPIXouter[5]={397,433,475,523,578};
3659
3660	if( r<12 ){
3661	if (j < MINPIXinner[r])
3662	{
3663	j = i+6 * r - int(rho_pixel);
3664	}
3665
3666	k=j-1;
3667	if (k < MINPIXinner[r])
3668	{
3669	k = MINPIXinner[r+1]-1;
3670	}
3671	}
3672	if( r > 11 && r < 16){
3673	if (j < MINPIXouter[r-12])
3674	{
3675	j = i + (r - 6) * 6 - int(rho_pixel);
3676
3677	}
3678
3679	k=j-1;
3680	if ( k < MINPIXouter[r-12])
3681	{
3682	k = MINPIXouter[r-11] - 1;
3683	}
3684	}
3685	rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k];
3686	}
3687	return (int(rho_pixel));
3688
3689	}
3690
3691
3692	//------------------------------------------------------------
3693	// @name produce_phes
3694	//
3695	// @desc read the photons of an event, pixelize them and simulate QE
3696	// @desc return various statistics and the array of Photoelectrons
3697	//
3698	// @date Mon Feb 14 16:44:21 CET 2000
3699	// @function @code
3700	//------------------------------------------------------------
3701
3702	int produce_phes( FILE *sp, // the input file
3703	class MGeomCam *camgeom, // the camera layout
3704	float minwl_nm, // the minimum accepted wavelength
3705	float maxwl_nm, // the maximum accepted wavelength
3706	class MTrigger *trigger, // the generated phes
3707	class MFadc *fadc,
3708	int *itotnphe, // total number of produced photoelectrons
3709	float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
3710	int *incph, // total number of cph read
3711	float *tmin_ns, // minimum arrival time of all phes
3712	float *tmax_ns, // maximum arrival time of all phes
3713	int ict // Telescope that is being analised to get the right QE.
3714	){
3715
3716	static uint i;
3717	static int k, ipixnum;
3718	static float cx, cy, wl, qe, t;
3719	static float cu, cv, cw;
3720	static MCCphoton photon;
3721	static float **qept;
3722	static char flag[SIZE_OF_FLAGS + 1];
3723	static float radius_mm;
3724
3725	// reset variables
3726
3727	for ( i=0; i<camgeom->GetNumPixels(); ++i ){
3728
3729	nphe[i] = 0.0;
3730
3731	}
3732
3733	*incph = 0;
3734
3735	radius_mm = camgeom->GetMaxRadius();
3736
3737
3738	//- - - - - - - - - - - - - - - - - - - - - - - - -
3739	// read photons and "map" them into the pixels
3740	//--------------------------------------------------
3741
3742	// initialize CPhoton
3743
3744	photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
3745
3746	// read the photons data
3747
3748
3749	// loop over the photons
3750
3751	while (1) {
3752
3753	fread ( flag, SIZE_OF_FLAGS, 1, sp );
3754
3755	if (isA( flag, FLAG_END_OF_EVENT ))
3756	break;
3757
3758	memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
3759
3760	fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
3761
3762
3763	// increase number of photons
3764
3765	(*incph)++;
3766
3767	// Chceck if photon is inside trigger time range
3768
3769	t = photon.get_t() ;
3770
3771	if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
3772	continue;
3773	}
3774
3775	//
3776	// Pixelization
3777	//
3778
3779	cx = photon.get_x();
3780	cy = photon.get_y();
3781
3782	// get wavelength
3783
3784	wl = photon.get_wl();
3785
3786	// cout << "wl " << wl << " radius " << sqrt(cxcx + cycy) << "\n";
3787
3788	// check if photon has valid wavelength and is inside outermost camera radius
3789
3790	if( (wl > maxwl_nm) \|\| (wl < minwl_nm) \|\| (sqrt(cxcx + cycy)*10 > radius_mm ) ){
3791	continue;
3792
3793	}
3794
3795	ipixnum = bpoint_is_in_pix(cx10, cy10, camgeom);
3796
3797	// -1 = the photon is in none of the pixels
3798	// 0 = the phton is in the central pixel, which is not used for trigger
3799	if (ipixnum==-1 \|\| ipixnum==0) {
3800	continue;
3801	}
3802
3803	//+++
3804	// QE simulation
3805	//---
3806
3807	// set pointer to the QE table of the relevant pixel
3808
3809	qept = (float **)QE[ict][ipixnum];
3810
3811	// check if wl is inside table; outside the table, QE is assumed to be zero
3812
3813	if((wl < qept[0][0]) \|\| (wl > qept[0][pointsQE[ict]-1])){
3814	continue;
3815
3816	}
3817
3818	// find data point in the QE table (-> k)
3819
3820	k = 1; // start at 1 because the condition was already tested for 0
3821	while (k < pointsQE[ict]-1 && qept[0][k] < wl){
3822	k++;
3823	}
3824
3825	// calculate the qe value between 0. and 1.
3826
3827	qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
3828
3829	// Apply incient angular correction due to Light Guides
3830
3831	cu=photon.get_u();
3832	cv=photon.get_v();
3833	cw=90.0-acos(sqrt(1-cucu-cvcv))*3.14159/180.0;
3834
3835	// find data point in the QE table (-> k)
3836
3837	k = 0;
3838	while (k < pointsWC-1 && WC[0][k] < cw){
3839	k++;
3840	}
3841
3842	// correct the qe with WC data.
3843
3844	qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
3845
3846	// if random > quantum efficiency, reject it
3847
3848	if ( (RandomNumber) > qe ) {
3849	continue;
3850	}
3851
3852	//+++
3853	// The photon has produced a photo electron
3854	//---
3855
3856	// cout << " accepted\n";
3857
3858	// increment the number of photoelectrons in the relevant pixel
3859
3860	nphe[ipixnum] += 1.0;
3861
3862	// store the new photoelectron
3863
3864	fadc->Fill(ipixnum,(t-*tmin_ns) ,
3865	trigger->FillShow(ipixnum,t-*tmin_ns),
3866	!((camgeom)[ipixnum].GetD()>(camgeom)[0].GetD()));
3867
3868	*itotnphe += 1;
3869	}
3870
3871	return(0);
3872
3873	}
3874
3875
3876	//------------------------------------------------------------
3877	// @name produce_nsbrates
3878	//
3879	// @desc read the starfield file, call produce_phes on it in,
3880	// @desc different wavebands, calculate the nsbrates
3881	//
3882	// @date Mon Feb 14 16:44:21 CET 2000
3883	// @function @code
3884	//------------------------------------------------------------
3885
3886	int produce_nsbrates( char *iname, // the starfield input file name
3887	MGeomCam *camgeom, // camera layout
3888	float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS], // the product of this function:
3889	// the NSB rates in phe/ns for each pixel
3890	int ict
3891	){
3892
3893	uint i, j; // counters
3894	int k, ii; // counters
3895
3896	static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
3897	WAVEBANDBOUND2,
3898	WAVEBANDBOUND3,
3899	WAVEBANDBOUND4,
3900	WAVEBANDBOUND5,
3901	WAVEBANDBOUND6 };
3902
3903	FILE *infile; // the input file
3904	fpos_t fileposition; // marker on the input file
3905	static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
3906	static MCEventHeader evth; // the event header
3907	static MCEventHeader evth_2; // the event header
3908	static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
3909	int reflector_file_version;
3910	int itnphe; // total number of produced photoelectrons
3911	int itotnphe; // total number of produced photoelectrons after averaging
3912	int incph; // total number of cph read
3913	float tmin_ns; // minimum arrival time of all phes
3914	float tmax_ns; // maximum arrival time of all phes
3915	float integtime_ns; // integration time from the starfield generator
3916	char flag_new[4];
3917
3918	// open input file
3919
3920	log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
3921
3922	infile = fopen( iname, "r" );
3923
3924	// check if the satrfield input file exists
3925
3926	if ( infile == NULL ) {
3927
3928	log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
3929
3930	log( SIGNATURE, "There is not NSB from the Stars\n");
3931
3932	return (0);
3933	}
3934
3935	// get signature, and check it
3936
3937	if((reflector_file_version=check_reflector_file( infile ))==FALSE){
3938	exit(1);
3939	}
3940
3941	// Instance of MTrigger and MFadc
3942
3943	MTrigger trigger(camgeom->GetNumPixels(),camgeom,Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
3944	MFadc flashadc(camgeom->GetNumPixels());
3945
3946	// initialize flag
3947
3948	strcpy( cflag, " \0" );
3949
3950	// get flag
3951
3952	fread( cflag, SIZE_OF_FLAGS, 1, infile );
3953
3954	if ( ! feof(infile)){
3955
3956	// reading .rfl file
3957
3958	if(!isA( cflag, FLAG_START_OF_RUN )){
3959	error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
3960	}
3961	else { // found start of run
3962
3963	if (reflector_file_version > 5){
3964
3965	fread( flag_new, 4, 1, infile );
3966
3967	if(!isA( flag_new, FLAG_START_OF_HEADER)){
3968
3969	// We break the main loop
3970	cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
3971	cout<<" We assume no Star Light"<<endl;
3972	return(0);
3973	}
3974
3975	Float_t flag_temp[SIZE_OF_MCRUNHEADER];
3976
3977	fread( flag_temp, (SIZE_OF_MCRUNHEADER)*sizeof(float), 1, infile );
3978
3979	}
3980
3981	fread( cflag, SIZE_OF_FLAGS, 1, infile );
3982
3983	if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
3984	fread( flag_new, 4, 1, infile );
3985
3986	if(!isA( flag_new, FLAG_EVENT_HEADER)){
3987
3988	// We break while events loop
3989	cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<endl;
3990	cout<<" We assume no light from Stars"<<endl;
3991	return(0);
3992	}
3993
3994	// get MCEventHeader
3995
3996	if (reflector_file_version<6)
3997	fread( (char*)&evth, evth.mysize(), 1, infile );
3998	else
3999	fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
4000
4001	if (reflector_file_version<6)
4002	integtime_ns = evth.get_energy();
4003	else
4004	integtime_ns = evth_2.get_energy();
4005
4006	// memorize where we are in the file
4007
4008	if (fgetpos( infile, &fileposition ) != 0){
4009	error( SIGNATURE, "Cannot position in file ...\n");
4010	}
4011
4012	// loop over the wavebands
4013
4014	for(i=0; i<iNUMWAVEBANDS; i++){
4015
4016	// initialize the rate array
4017
4018	for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4019	rate_phepns[j][i] = 0.;
4020	}
4021
4022	itotnphe = 0;
4023
4024	// read the photons and produce the photoelectrons
4025	// - in order to average over the QE simulation, call the
4026	// production function iNUMNSBPRODCALLS times
4027
4028	for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
4029
4030	// position the file pointer to the beginning of the photons
4031
4032	fsetpos( infile, &fileposition);
4033
4034	// produce photoelectrons
4035
4036	k = produce_phes( infile,
4037	camgeom,
4038	wl_nm[i],
4039	wl_nm[i+1],
4040	&trigger, // this is a dummy here
4041	&flashadc, // this is a dummy here
4042	&itnphe,
4043	nphe, // we want this!
4044	&incph,
4045	&tmin_ns,
4046	&tmax_ns,
4047	0
4048	); // photons from starfield
4049
4050	if( k != 0 ){ // non-zero returnvalue means error
4051	cout << "Exiting.\n";
4052	exit(1);
4053	}
4054
4055	for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4056	rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
4057	}
4058
4059	itotnphe += itnphe;
4060
4061	} // end for(ii=0 ...
4062
4063	fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
4064	wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
4065
4066	} // end for(i=0 ...
4067
4068	}
4069	else{
4070	cout << "Starfield file contains no event.\nExiting.\n";
4071	exit(1);
4072	} // end if( isA ... event
4073	} // end if ( isA ... run
4074	}
4075	else{
4076	cout << "Starfield file contains no run.\nExiting.\n";
4077	exit(1);
4078	}
4079
4080	fclose( infile );
4081
4082	return(0);
4083
4084	}
4085
4086
4087	//------------------------------------------------------------
4088	// @name produce_nsb_phes
4089	//
4090	// @desc produce the photoelectrons from the nsbrates
4091	//
4092	// @date Thu Feb 17 17:10:40 CET 2000
4093	// @function @code
4094	//------------------------------------------------------------
4095
4096	int produce_nsb_phes( float *atmin_ns,
4097	float *atmax_ns,
4098	float theta_rad,
4099	struct camera *cam,
4100	float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
4101	float difnsbr_phepns[iMAXNUMPIX],
4102	float extinction[iNUMWAVEBANDS],
4103	float fnpx[iMAXNUMPIX],
4104	MTrigger *trigger,
4105	MFadc *fadc,
4106	int *inphe,
4107	float base_mv[iMAXNUMPIX]){
4108
4109	float simtime_ns; // NSB simulation time
4110	int i, j, k, ii;
4111	float zenfactor; // correction factor calculated from the extinction
4112	int inumnsbphe; // number of photoelectrons caused by NSB
4113
4114	float t;
4115	TRandom random; // Random numbers generator
4116
4117	random.SetSeed(get_seeds(1));
4118
4119	ii = *inphe; // avoid dereferencing
4120
4121	// check if the arrival times are set; if not generate them
4122
4123	if(atmin_ns <SIMTIMEOFFSET_NS \|\| atmin_ns > *atmax_ns){
4124	*atmin_ns = 0.;
4125	*atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
4126
4127	}
4128	else{ // extend the event time window by the given offsets
4129
4130	atmin_ns = atmin_ns - SIMTIMEOFFSET_NS;
4131	atmax_ns = atmax_ns + SIMTIMEOFFSET_NS;
4132
4133	simtime_ns = atmax_ns - atmin_ns;
4134
4135	// make sure the simulated time is long enough for the FADC
4136	// simulation and not too long
4137
4138	if(simtime_ns< TOTAL_TRIGGER_TIME){
4139	atmin_ns = atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
4140	atmax_ns = atmin_ns + TOTAL_TRIGGER_TIME;
4141	simtime_ns = TOTAL_TRIGGER_TIME;
4142	}
4143
4144	if(simtime_ns> TOTAL_TRIGGER_TIME){
4145	atmax_ns =atmin_ns + TOTAL_TRIGGER_TIME;
4146	simtime_ns = TOTAL_TRIGGER_TIME;
4147	}
4148
4149	}
4150
4151	// initialize baselines
4152
4153	for(i=0; i<cam->inumpixels; i++){
4154	base_mv[i] = 0.;
4155	}
4156
4157	// calculate baselines and generate phes
4158
4159	for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
4160
4161	// calculate the effect of the atmospheric extinction
4162
4163	zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
4164
4165	for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
4166
4167	inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
4168	* simtime_ns );
4169
4170	base_mv[j] += inumnsbphe;
4171
4172	// randomize
4173
4174	if (inumnsbphe>0.0){
4175	inumnsbphe = random.Poisson (inumnsbphe );
4176	}
4177
4178	// create the photoelectrons
4179
4180	for(k=0; k < inumnsbphe; k++){
4181
4182	t=(RandomNumber * simtime_ns);
4183
4184	(fadc).Fill(j,t ,(trigger).FillNSB(j,t));
4185
4186	ii++; // increment total number of photoelectons
4187
4188	fnpx[j] += 1.; // increment number of photoelectrons in this pixel
4189
4190	}
4191
4192	} // end for(j=0 ...
4193	} // end for(i=0 ...
4194
4195	// finish baseline calculation
4196
4197	for(i=0; i<cam->inumpixels; i++){
4198	base_mv[i] = RESPONSE_FWHM RESPONSE_AMPLITUDE / simtime_ns;
4199	}
4200
4201	*inphe = ii; // update the pointer
4202
4203	return(0);
4204	}
4205
4206
4207	// @endcode
4208
4209
4210	//=------------------------------------------------------------
4211	//!@subsection Log of this file.
4212
4213	//!@{
4214	//
4215	// $Log: not supported by cvs2svn $
4216	// Revision 1.65 2003/10/26 19:43:00 blanch
4217	// - The screen output information has been improved to prevent some
4218	// non-desired running conditions just looking at the output screen.
4219	// - One MMcEvt for each Telscope is stored in the output file.
4220	// - 500 empty events are simualted to get a more precise estimation of the
4221	// pedestal Sigma for each pixel.
4222	//
4223	// Revision 1.64 2003/10/21 07:42:50 blanch
4224	// A factor 2.35 to transform the fwhm into the sigma of gaussian was missing
4225	// in the storing of FADC single hpe pulse determination.
4226	//
4227	// Revision 1.63 2003/10/17 19:38:31 blanch
4228	// Now the camera program will stop if a undefined Geometry is required.
4229	// The NSB is internally scaled for any camera geometry and qe.
4230	// The seeds in the input card are used to initilize the random numbers.
4231	// The Amplitud stored in the MMcFadcHeader is the amplitud of the sphe reponse.
4232	// The Pedestal rms is simualted in a artificail empty event.
4233	//
4234	// Revision 1.62 2003/09/26 11:25:07 blanch
4235	// Modification to be able to read MGeomCam branch for any Geometry.
4236	//
4237	// Revision 1.61 2003/09/25 17:09:20 blanch
4238	// Bug on the number of phe from diffuse NSB fixed.
4239	//
4240	// Revision 1.60 2003/09/23 16:50:55 blanch
4241	// WE do not read ct_file anymore since all Telescope information is
4242	// in the reflector or in MGeomCam.
4243	//
4244	// Revision 1.58 2003/09/15 09:59:53 blanch
4245	// The concept of the camera prgoram has not changed but this version has
4246	// quite a lot of changes to allow several Camera geometries as well as
4247	// multitelescope.
4248	//
4249	// It is suposed to be the first working code for camera 0.7.
4250	//
4251	// Revision 1.57 2003/07/17 18:02:46 blanch
4252	// Several new features introduced as well as fixed bugs
4253	//
4254	// - 1/100 events printed out
4255	// - Low gain implemented
4256	// - Different response for outer and inner pixels
4257	// - Some warnings removed
4258	// - pedestal and qe file from inpuit card
4259	// - Faster electronic simulation
4260	//
4261	// Revision 1.55 2003/02/12 12:22:10 blanch
4262	// * empty log message *
4263	//
4264	// Revision 1.54 2003/02/12 11:55:01 blanch
4265	// * empty log message *
4266	//
4267	// Revision 1.53 2003/01/23 18:35:21 blanch
4268	// * empty log message *
4269	//
4270	// Revision 1.52 2003/01/20 17:19:20 blanch
4271	// It fills the MMcCorsikaRun.
4272	//
4273	// Revision 1.51 2003/01/14 13:40:17 blanch
4274	// MRawRunHeader::fNumEvents has been filled with the number of events in
4275	// this file.
4276	// Problems in fImpact computation have been solved.
4277	// Option to set a dc value to rise the discriminator threshold has been added.
4278	//
4279	// Revision 1.50 2003/01/07 16:33:31 blanch
4280	// Star Field Rotation has been implemented by Raul Orduna. Now there is a
4281	// rotation for each shower. It is done by a non enter pixel rotation assuming
4282	// a circular symetry of the camera. It is not exact but it is accurate enough and
4283	// much faster.
4284	//
4285	// Revision 1.49 2002/12/13 10:04:07 blanch
4286	// * empty log message *
4287	//
4288	// Revision 1.48 2002/12/12 17:40:50 blanch
4289	// * empty log message *
4290	//
4291	// Revision 1.47 2002/12/10 17:19:31 blanch
4292	// * empty log message *
4293	//
4294	// Revision 1.46 2002/11/08 17:51:00 blanch
4295	// * empty log message *
4296	//
4297	// Revision 1.45 2002/10/29 17:15:27 blanch
4298	// Reading several reflector versions.
4299	//
4300	// Revision 1.44 2002/10/18 16:53:03 blanch
4301	// Modification to read several reflector files.
4302	//
4303	// Revision 1.43 2002/09/13 10:53:39 blanch
4304	// Minor change to remove some undisired comments.
4305	//
4306	// Revision 1.42 2002/09/09 16:00:49 blanch
4307	// Statement has been included to avoid writting to disk MParContainer and MArray.
4308	// It has also been added the effect of the WC, the actual values must be added,
4309	// once they are measured.
4310	//
4311	// Revision 1.41 2002/09/04 09:57:42 blanch
4312	// Modifications done to use MGeomCam from MARS.
4313	//
4314	// Revision 1.40 2002/07/16 16:15:22 blanch
4315	// A first implementation of the Star field rotation has been done.
4316	//
4317	// Revision 1.39 2002/04/27 10:48:39 blanch
4318	// Some unused varibles have been removed.
4319	//
4320	// Revision 1.38 2002/03/18 18:44:29 blanch
4321	// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
4322	//
4323	// Revision 1.37 2002/03/18 16:42:20 blanch
4324	// The data member fProductionSite of the MMcRunHeader has been set to 0,
4325	// which means that the prodution site is unknown.
4326	//
4327	// Revision 1.35 2002/03/15 15:15:52 blanch
4328	// Several modifications to simulate the actual trigger zone.
4329	//
4330	// Revision 1.34 2002/03/13 18:13:56 blanch
4331	// Some changes to fill correctly the new format of MMcRunHeader.
4332	//
4333	// Revision 1.33 2002/03/04 17:21:48 blanch
4334	// Small and not important changes.
4335	//
4336	// Revision 1.32 2002/02/28 15:04:52 blanch
4337	// A small back has been solved. Before, while not using the option
4338	// writte_all_images, not all triggered showers were stored. Now it is solved.
4339	// For that it is importatn taht the less restrictive trigger option is
4340	// checked first.
4341	// A new facility has been introduced and now one can choose the step size in
4342	// trigger loop mode for the discriminator threshold.
4343	// The close-compact topology for less than 3 pixels does not make sense. Before
4344	// the program was ignoring that, now it switch to simple neighbour condition.
4345	//
4346	// Revision 1.31 2002/01/18 17:41:02 blanch
4347	// The option of adding noise to all pixels or to not adding the noise
4348	// has been added.
4349	// We removed the pixels larger than 577. When there were more than one
4350	// trigger in one shower, the pixel number was increasing. Now it is
4351	// flagged by the variable MMcTrig::fFirstLvlTrig.
4352	//
4353	// Revision 1.30 2001/11/27 09:49:54 blanch
4354	// Fixing bug which was treating wrongly the extension of star photons.
4355	//
4356	// Revision 1.29 2001/11/14 17:38:23 blanch
4357	// Sveral changes have been done:
4358	// - bpoint_in_in_pixel has been dodified to speed up the program
4359	// - Some minor changes have been done to remove warnings
4360	// - buffer size and split version of the Branches have been removed
4361	// - Some modifications were needed to adat the program to the new
4362	// MRawEvtData::DeletePixels
4363	//
4364	// Revision 1.28 2001/10/26 16:31:45 blanch
4365	// Removing several warnings.
4366	//
4367	// Revision 1.27 2001/09/05 10:04:33 blanch
4368	// * empty log message *
4369	//
4370	// Revision 1.26 2001/07/19 09:29:53 blanch
4371	// Different threshold for each pixel can be used.
4372	//
4373	// Revision 1.25 2001/05/08 08:07:54 blanch
4374	// New numbering for branches from different trigger conditions has been
4375	// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
4376	// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
4377	// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
4378	// information already filled.
4379	//
4380	// Revision 1.24 2001/04/06 16:48:09 magicsol
4381	// New camera version able to read the new format of the reflector output:
4382	// reflector 0.4
4383	//
4384	// Revision 1.23 2001/03/28 16:13:41 blanch
4385	// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
4386	// been solved.
4387	//
4388	// Revision 1.22 2001/03/20 18:52:43 blanch
4389	// * empty log message *
4390	//
4391	// Revision 1.21 2001/03/19 19:53:03 blanch
4392	// Some print outs have been removed.
4393	//
4394	// Revision 1.20 2001/03/19 19:30:06 magicsol
4395	// Minor changes have been done to improve the FADC pedestals treatment.
4396	//
4397	// Revision 1.19 2001/03/05 11:14:41 magicsol
4398	// I changed the position of readinf a parameter. It is a minnor change.
4399	//
4400	// Revision 1.18 2001/03/05 10:36:52 blanch
4401	// A branch with information about the FADC simulation (MMcFadcHeader) is writen
4402	// in the McHeaders tree of the aoutput root file.
4403	// The NSB contribution is only applied if the the number of phe form the shower
4404	// are above a given number.
4405	//
4406	// Revision 1.17 2001/02/23 11:05:57 magicsol
4407	// Small changes due to slightly different output format and the introduction of
4408	// pedesals for teh FADC.
4409	//
4410	// Revision 1.16 2001/01/18 18:44:40 magicsol
4411	// * empty log message *
4412	//
4413	// Revision 1.15 2001/01/17 09:32:27 harald
4414	// The changes are neccessary to have the same name for trees in MC and in
4415	// data. So now there should be no differences in MC data and real data from
4416	// FADC system.
4417	//
4418	// Revision 1.14 2001/01/15 12:33:34 magicsol
4419	// Some modifications have been done to use the new (Dec'2000) Raw data format.
4420	// There are also some minnors modifications to adapt the improvements in the
4421	// MTrigger class (overlaping time and trigger cells).
4422	//
4423	// Revision 1.13 2000/10/25 08:14:23 magicsol
4424	// The routine that produce poisson random numbers to decide how many phe
4425	// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
4426	// is used.
4427	//
4428	// Revision 1.10 2000/07/04 14:10:20 MagicSol
4429	// Some changes have been done in the root output file. The RawEvt tree is only
4430	// stored in the single trigger mode.
4431	// The trigger input parameters are also given by the general input card.
4432	// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
4433	// each one can be studied seperately.
4434	//
4435	// Revision 1.9 2000/06/13 13:25:24 blanch
4436	// The multiple arrays have been replaced, since they do not work
4437	// in alpha machines. Now we are using pointers and the command new
4438	// to allocate memory.
4439	//
4440	// Revision 1.8 2000/05/11 13:57:27 blanch
4441	// The option to loop over several trigger configurations has been included.
4442	// Some non-sense with triggertime range has been solved.
4443	// Montecarlo information and ADC counts are saved in a root file.
4444	// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
4445	//
4446	// Revision 1.7 2000/03/24 18:10:46 blanch
4447	// A first FADC simulation and a trigger simulation are already implemented.
4448	// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
4449	// 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.
4450	//
4451	// Revision 1.6 2000/03/20 18:35:11 blanch
4452	// 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.
4453	//
4454	// Revision 1.5 2000/02/18 17:40:35 petry
4455	// This version includes drastic changes compared to camera.cxx 1.4.
4456	// It is not yet finished and not immediately useful because the
4457	// trigger simulation is not yet re-implemented. I had to take it
4458	// out together with some other stuff in order to tidy the whole
4459	// program up. This is not meant as an insult to anyone. I needed
4460	// to do this in order to be able to work on it.
4461	//
4462	// This version has been put in the repository in order to be
4463	// able to share the further development with others.
4464	//
4465	// If you need something working, wait or take an earlier one.
4466	// See file README.
4467	//
4468	// Revision 1.4 2000/01/25 08:36:23 petry
4469	// The pixelization in previous versions was buggy.
4470	// This is the first version with a correct pixelization.
4471	//
4472	// Revision 1.3 2000/01/20 18:22:17 petry
4473	// Found little bug which makes camera crash if it finds a photon
4474	// of invalid wavelength. This bug is now fixed and the range
4475	// of valid wavelengths extended to 290 - 800 nm.
4476	// This is in preparation for the NSB simulation to come.
4477	// Dirk
4478	//
4479	// Revision 1.2 1999/11/19 08:40:42 harald
4480	// Now it is possible to compile the camera programm under osf1.
4481	//
4482	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
4483	// This the starting point for CVS controlled further developments of the
4484	// camera program. The program was originally written by Jose Carlos.
4485	// But here you can find a "rootified" version to the program. This means
4486	// that there is no hbook stuff in it now. Also the output of the
4487	// program changed to the MagicRawDataFormat.
4488	//
4489	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
4490	//
4491	// Revision 1.3 1999/10/22 15:01:28 petry
4492	// version sent to H.K. and N.M. on Fri Oct 22 1999
4493	//
4494	// Revision 1.2 1999/10/22 09:44:23 petry
4495	// first synthesized version which compiles and runs without crashing;
4496	//
4497	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
4498	// first synthesised version
4499	//
4500	// Revision 1.13 1999/03/15 14:59:05 gonzalez
4501	// camera-1_1
4502	//
4503	// Revision 1.12 1999/03/02 09:56:10 gonzalez
4504	// * empty log message *
4505	//
4506	//
4507	//!@}
4508
4509	//=EOF

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

Download in other formats: