Context Navigation

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

Visit:

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

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

Download in other formats: