Context Navigation

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

Visit:

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

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

Download in other formats: