Context Navigation

camera.cxx@ 2964

Visit:

Last change on this file since 2964 was 2962, checked in by blanch, 21 years ago

[Changes mainly done by A. Moralejo] Several minnor changes have been done. For instance, some name of the variables have been modified to a more self-explained name and modifications while reading the asciis files at the end of the reflector file has been introduced. The possibilty to unlarge the point spread function has been introduced in order to be able to simualte the current data. All pixels are always written.

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

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

Context Navigation

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

Download in other formats: