Context Navigation

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

Visit:

Last change on this file since 6569 was 6567, checked in by moralejo, 21 years ago
Corrected bug in the setting of the trigger patterns. It turns out that in the current data format, MRawRunHeader.fFormatVersion=5, as of February 2005, the trigger bits are "inverted" ( 0 <-> 1 ) in their meaning (see Mars/mtrigger/MTriggerPattern).
File size: 163.4 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	//!@section Source code of \|camera.cxx\|.
16
17	//=-----------------------------------------------------------
18	//!@subsection Includes and Global variables definition.
19
20	//!@{
21
22	// includes for ROOT
23	// BEWARE: the order matters!
24
25	#include "TROOT.h"
26
27	#include "TRandom.h"
28	#include "TApplication.h"
29
30	#include "TFile.h"
31	#include "TTree.h"
32	#include "TBranch.h"
33	#include "TCanvas.h"
34
35	#include "TArrayC.h"
36	#include "TObjArray.h"
37
38	#include "MTrigger.hxx"
39	#include "MFadc.hxx"
40	#include "MLons.hxx"
41
42	#include "MRawRunHeader.h"
43	#include "MRawEvtData.h"
44	#include "MRawEvtHeader.h"
45	#include "MRawCrateArray.h"
46	#include "MRawCrateData.h"
47	#include "MMcEvt.hxx"
48	#include "MMcEvtBasic.h"
49	#include "MMcTrig.hxx"
50	#include "MMcRunHeader.hxx"
51	#include "MMcConfigRunHeader.h"
52	#include "MMcCorsikaRunHeader.h"
53	#include "MMcTrigHeader.hxx"
54	#include "MMcFadcHeader.hxx"
55	#include "MGeomCamMagic.h"
56	#include "MGeomCamMagic919.h"
57	#include "MGeomCamMagicHG.h"
58	#include "MGeomCamECO1000.h"
59	#include "MGeomCamECO1000HG.h"
60	#include "MGeomCamCT1.h"
61	#include "MGeomCamCT1Daniel.h"
62	#include "MGeomPix.h"
63	#include "MGeomCorsikaCT.h"
64	#include "MTriggerPattern.h"
65
66	/*!@"
67
68	All the defines are located in the file \|camera.h\|.
69
70	@"*/
71
72	#include "camera.h"
73
74	//!@}
75
76	/*!@"
77
78	The following set of flags are used in time of compilation. They do
79	not affect directly the behaviour of the program at run-time
80	(though, of course, if you disconnected the option for
81	implementation of the Trigger logic, you will not be able to use any
82	trigger at all. The 'default' values mean default in the sense of
83	what you got from the server when you obtained this program.
84
85	@"*/
86
87	//!@{
88
89	// flag for debugging (default: OFF )
90	#define __DEBUG__
91	#undef __DEBUG__
92
93
94	//!@}
95
96	//=-----------------------------------------------------------
97	//!@subsection Definition of global variables.
98
99	/*!@"
100
101	Now we define some global variables with data about the telescope,
102	such as "focal distance", number of pixels/mirrors,
103	"size of the camera", and so on.
104
105	@"*/
106
107	/*!@"
108
109	Depending on the telescope we are using (CT1 or MAGIC), the
110	information stored in the definition file is different.
111	The variable \|ct_Type\| has the value 0 when we use
112	CT1, and 1 when we use MAGIC.
113
114	@"*/
115
116	/*!@"
117
118	And this is the information about the whole telescope.
119
120	@"*/
121
122	//!@{
123
124	// parameters of the CT (from the CT definition file)
125
126	//@: Number of pixels
127	static int ct_NPixels;
128
129	//@: Number of CT
130	static int ct_Number;
131
132	//@: list of showers to be skipped
133	static int *Skip;
134
135	//@: number of showers to be skipped
136	static int nSkip=0;
137
138	//@: flag: TRUE: data come from STDIN; FALSE: from file
139	static int Data_From_STDIN = FALSE;
140
141	//@: flag: TRUE: write the event headers of all events to output;
142	//@: FALSE: only triggered showers
143	static int Write_All_Event_Headers = TRUE;
144
145	static int Write_McEvt = TRUE;
146	static int Write_McTrig = TRUE;
147	static int Write_McFADC = TRUE;
148	static int Write_RawEvt = FALSE;
149
150	//@: flag: TRUE: selection on the energy
151	static int Select_Energy = TRUE;
152
153	//@: Lower edge of the selected energy range (in GeV)
154	static float Select_Energy_le = 0.0;
155
156	//@: Upper edge of the selected energy range (in GeV)
157	static float Select_Energy_ue = 100000.0;
158
159	//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
160	static int FADC_Scan = FALSE;
161
162	//@: flag: TRUE: show all trigger signal in the screen; FALSE: don't
163	static int Trigger_Scan = FALSE;
164
165	//@: flag: TRUE: loop trigger analysis over several thresholds, multiplicities and topologies; FALSE: a single trigger configuration
166	static int Trigger_Loop = FALSE;
167
168	//@: flag: TRUE: Different threshold for each pixel ; FALSE: same threshold for all pixels
169	static int Individual_Thres_Pixel = FALSE;
170
171	//@: Properties of the trigger
172	static float Trigger_gate_length = 6.0;
173	static float Trigger_response_ampl = 1.0;
174	static float Trigger_response_fwhm = 2.0;
175	static float Trigger_overlaping_time= 0.25;
176	static float Trigger_noise= 0.3;
177
178	//@: Properties of the FADC
179	static Int_t FADC_shape = 0;
180	static float FADC_response_integ = MFADC_RESPONSE_INTEGRAL;
181	static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
182	static Int_t FADC_shape_out = 0;
183	static float FADC_resp_integ_out = MFADC_RESPONSE_INTEGRAL;
184	static float FADC_resp_fwhm_out = MFADC_RESPONSE_FWHM;
185	static float FADC_slices_per_ns = FADC_SLICES_PER_NSEC;
186	static Int_t FADC_slices_written = FADC_SLICES;
187	static float FADC_noise_inner = 2.0;
188	static float FADC_noise_outer = 2.0;
189	static float DIGITAL_noise = 0.0;
190	static float FADC_high2low = HIGH2LOWGAIN;
191
192	//@: Trigger conditions for a single trigger mode
193	static float **qThreshold;
194	static int Trigger_multiplicity[MAX_NUMBER_OF_CTS];
195	static int Trigger_topology[MAX_NUMBER_OF_CTS];
196
197	//@: Upper and lower edges of the trigger loop
198	static float Trigger_loop_lthres = 2.0;
199	static float Trigger_loop_uthres = 10.0;
200	static float Trigger_loop_sthres = 1.0;
201	static int Trigger_loop_lmult = 2;
202	static int Trigger_loop_umult = 10;
203	static int Trigger_loop_ltop = 0;
204	static int Trigger_loop_utop = 2;
205
206	//@: Direction of each shower
207	static float Zenith = 0.0;
208	static float Azimutal = 90.0;
209
210	//@: Mispointing Simulation
211	static int missPointing = 0;
212	static float missP_x = 0.0;
213	static float missP_y = 0.0;
214
215	//@: Point Spread Function Added
216	static float Spot_x=0.0;
217	static float Spot_y=0.0;
218	static float Spotsigma=0.0;
219
220	//@: PMT time jitter
221	static float pmt_jitter=0.25;
222
223
224	//!@}
225
226	/*!@"
227
228	The following double-pointer is a 2-dimensional table with information
229	about each pixel. The routine read_pixels will generate
230	the information for filling it using igen_pixel_coordinates().
231
232	@"*/
233
234	//!@{
235	// Pointer to a tables/Arrays with information about the pixels
236	// and data stored on them with information about the pixels
237
238
239	//@: contents of the pixels (ph.e.)
240	static float *fnpix;
241
242
243	//!@}
244
245	/*!@"
246
247	The following double-pointer is a 2-dimensional table with the
248	Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
249	of the wavelength @$\lambda@$. The routine \|read_pixels()\| will read
250	also this information from the file \|qe.dat\|.
251
252	@"*/
253
254	//!@{
255	// Pointer to a table with QE, number of datapoints, and wavelengths
256
257	//@: table of QE
258	static float ****QE;
259
260	//@: number of datapoints for the QE curve
261	static int pointsQE[MAX_NUMBER_OF_CTS];
262
263	//@: table of QE
264	static float *QElambda;
265
266	//@: table of lightguide = WC; WC_outer for outer pixels.
267	static float **WC;
268	static float **WC_outer;
269
270	//@: number of datapoints for the WC curve
271	static int pointsWC;
272
273	//!@}
274
275	/*!@"
276
277	The following double-pointer is a 2-dimensional table with information
278	about each mirror in the dish. The routine \|read_ct_file()\| will read
279	this information from the CT definition file.
280
281	@"*/
282
283	/*!@"
284
285	We define a table into where random numbers will be stored.
286	The routines used for random number generation are provided by
287	\|RANLIB\| (taken from NETLIB, \|www.netlib.org\|), and by
288	the routine \|double drand48(void)\| (prototype defined in
289	\|stdlib.h\|) through the macro \|RandomNumber\| defined in
290	\|camera.h\|.
291
292	@"*/
293
294
295	//!@}
296
297	extern char FileName[];
298	// switch on starfield rotation
299	static int Starfield_rotate = TRUE;
300
301	//=-----------------------------------------------------------
302	// @subsection Main program.
303
304	//!@{
305
306	//++++++++++++++++++++++++++++++++++++++++
307	// MAIN PROGRAM
308	//----------------------------------------
309
310	int main(int argc, char **argv)
311	{
312
313	//!@' @#### Definition of variables.
314	//@'
315
316	char **inname_CT; //@< array of names for each CT input file
317	char starfieldname[256]; //@< starfield input file name
318	char qe_filename[256]; //@< qe input file name
319
320	char datname[256]; //@< data (ASCII) output file name
321
322	char rootname[256] ; //@< ROOT file name
323	char rootname_loop[256] ; //@< ROOT file name
324
325	char parname[256]; //@< parameters file name
326
327	char nsbpathname[256]; //@< directory with the dataabse for th NSB
328	char nsbpath_outer[256]; //@< directory with the dataabse for outer pixels
329
330	char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
331	char flag_new[SIZE_OF_FLAGS + 1]; //@< New flag for the run header in the .rfl file
332	char **GeometryName; //@< Name of MGeomCam for each CT
333	int GeometryCamera[MAX_NUMBER_OF_CTS]; //@< Identification of MGeomCam for each CT
334	int TriggerPixels[MAX_NUMBER_OF_CTS]; //@< Number of pixels in the trigger region for each CT
335
336	int reflector_file_version=0; //@< vrsion of he reflector file
337
338	FILE *inputfile[MAX_NUMBER_OF_CTS]; //@< stream for the input file
339
340	ofstream datafile; //@< stream for the data file
341
342	MCRunHeader mcrunh; //@< Run Header class (MC)
343	MCEventHeader mcevth[MAX_NUMBER_OF_CTS]; //@< Event Header class (MC)
344	MCEventHeader_2 mcevth_2[MAX_NUMBER_OF_CTS]; //@< Event Header class (MC) for reflector > 0.6
345	MCCphoton cphoton; //@< Cherenkov Photon class (MC)
346
347	int inumphe; //@< number of photoelectrons in an event from showers
348	int inumphe_CT[MAX_NUMBER_OF_CTS]; //@< number of photoelectrons in an event from showers
349	float inumphensb[MAX_NUMBER_OF_CTS]; //@< number of photoelectrons in an event from nsb
350
351	float arrtmin_ns; //@ arrival time of the first photoelectron
352	float arrtmax_ns; //@ arrival time of the last photoelectron
353
354	float thetaCT[MAX_NUMBER_OF_CTS], phiCT[MAX_NUMBER_OF_CTS]; //@< theta and phi of telescopes
355
356	//@: Coordinates of telescopes in Corsika's coordinate system
357	float CTx[MAX_NUMBER_OF_CTS];
358	float CTy[MAX_NUMBER_OF_CTS];
359	float CTz[MAX_NUMBER_OF_CTS];
360
361	float mirror_frac[MAX_NUMBER_OF_CTS];
362
363	float thetashw, phishw; //@< parameters of a given shower
364	float coreX, coreY; //@< core position
365	float impactD[MAX_NUMBER_OF_CTS]; //@< impact parameter
366	float l1, m1, n1; //@< auxiliary variables
367	float factorqe_NSB[MAX_NUMBER_OF_CTS]; //@< factor on the NSB depending of QE
368
369	int nshow=0; //@< partial number of shower in a given run
370	int ntshow=0; //@< total number of showers
371	int ncph[MAX_NUMBER_OF_CTS]; //@< partial number of photons in a given run
372	int ntcph[MAX_NUMBER_OF_CTS]; //@< total number of photons
373
374	int ibr, j, k; //@< simple counters
375
376	int addElecNoise; //@< Will we add ElecNoise?
377
378	float riseDiskThres;
379	float secureDiskThres;
380
381	int simulateNSB; //@< Will we simulate NSB?
382	int nphe2NSB; //@< From how many phe will we simulate NSB?
383	float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel)
384	float **diffnsb_phepns;
385	//@< diffuse NSB values for each pixel derived from meanNSB
386
387	float **nsbrate_phepns; //@< non-diffuse nsb
388	//@< photoelectron rates
389	float **nsb_phepns;
390	//@< NSB values for each pixel
391
392	float **nsb_phepns_rotated;
393	//@< NSB values for each pixel after rotation
394
395	Float_t nsb_trigresp[TRIGGER_TIME_SLICES]; //@< array to write the trigger
396	//@< response from the database
397	Float_t *nsb_fadcresp; //@< array to write the fadc
398	//@< response from the database
399	Byte_t trigger_map[((Int_t)(CAMERA_PIXELS/8))+1]; //@< Pixels on when the
400	//@< camera triggers
401	Float_t fadc_elecnoise[CAMERA_PIXELS]; //@< Electronic noise for each pixel
402	Float_t fadc_diginoise[CAMERA_PIXELS]; //@< Digital noise for each pixel
403	Float_t fadc_pedestals[CAMERA_PIXELS]; //@< array for fadc pedestals values
404	Float_t fadc_sigma[CAMERA_PIXELS]; //@< array for fadc pedestals sigma
405	Float_t fadc_sigma_low[CAMERA_PIXELS]; //@< array for fadc pedestals sigma
406
407	float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
408	EXTWAVEBAND1,
409	EXTWAVEBAND2,
410	EXTWAVEBAND3,
411	EXTWAVEBAND4,
412	EXTWAVEBAND5
413	};
414	float zenfactor=1.0; // correction factor calculated from the extinction
415
416	int nstoredevents = 0;
417	int flagstoring = 0;
418
419	int ntrigger[MAX_NUMBER_OF_CTS]; //@< number of triggers in the whole file
420	int btrigger = 0; //@< trigger flag
421	int ithrescount; //@< counter for loop over threshold trigger
422	float fthrescount; //@< value for loop over threshold trigger
423	int imulticount; //@< counter for loop over multiplicity trigger
424	int itopocount; //@< counter for loop over topology trigger
425	int isorttopo[3]; //@< sorting the topologies
426	int icontrigger; //@< number of trigger conditions to be analised
427
428	float fpixelthres[CAMERA_PIXELS]; //@< Threshold values
429
430	TArrayC *fadcValues; //@< the analog Fadc High gain signal for pixels
431	TArrayC *fadcValuesLow; //@< the analog Fadc Low gain signal for pixels
432
433	int still_in_loop = FALSE;
434
435	//@< Variables to fill the McRunHeader
436	Int_t sfRaH = 5;
437	Int_t sfRaM = 34;
438	Int_t sfRaS = 32;
439	Int_t sfDeD = 22;
440	Int_t sfDeM = 00;
441	Int_t sfDeS = 52;
442	Float_t shthetamax = 0.0;
443	Float_t shthetamin = 0.0;
444	Float_t shphimax = 0.0;
445	Float_t shphimin = 0.0;
446	UInt_t corsika = 5200 ;
447	Float_t maxpimpact = 0.0;
448
449	// Star Field Rotation variables
450	Float_t zenith = 0.0;
451	Float_t azimutal = 90.0;
452	Float_t rho , C3 , C2 , C1;
453
454	//!@' @#### Definition of variables for \|getopt()\|.
455	//@'
456
457	int ch, errflg = 0; //@< used by getopt
458
459	//!@' @#### Initialisation of some variables
460	//@'
461
462	inname_CT = new char *[MAX_NUMBER_OF_CTS];
463	GeometryName = new char *[MAX_NUMBER_OF_CTS];
464	diffnsb_phepns = new float *[MAX_NUMBER_OF_CTS];
465	nsb_phepns = new float *[MAX_NUMBER_OF_CTS];
466	nsb_phepns_rotated = new float *[MAX_NUMBER_OF_CTS];
467
468	for(int i=0;i<MAX_NUMBER_OF_CTS;i++)
469	{
470	inname_CT[i] = new char[256];
471	GeometryName[i] = new char[50];
472	diffnsb_phepns[i] = new float[iMAXNUMPIX];
473	nsb_phepns[i] = new float[iMAXNUMPIX];
474	nsb_phepns_rotated[i] = new float[iMAXNUMPIX];
475	CTx[i] = 0.; CTy[i] = 0.; CTz[i] = 0.;
476	}
477
478	nsbrate_phepns = new float *[iMAXNUMPIX];
479	for(int i=0;i<iMAXNUMPIX;i++)
480	nsbrate_phepns[i] = new float[iNUMWAVEBANDS];
481
482	qThreshold = new float *[MAX_NUMBER_OF_CTS];
483	for(int i=0;i<MAX_NUMBER_OF_CTS;i++)
484	qThreshold[i] = new float [CAMERA_PIXELS];
485
486	for(int i=0;i<iMAXNUMPIX;i++){
487	for(int ict=0;ict<MAX_NUMBER_OF_CTS;ict++){
488	nsb_phepns[ict][i]=0;
489	nsb_phepns_rotated[ict][i]=0;
490	}
491	for(j=0;j<iNUMWAVEBANDS;j++)
492	nsbrate_phepns[i][j]=0.0; //@< Starlight rates initialised at 0.0
493	}
494	for(int i=0;i<MAX_NUMBER_OF_CTS;i++)
495	ntcph[i]=0;
496	/*!@'
497
498	@#### Beginning of the program.
499
500	We start with the main program. First we (could) make some
501	presentation, and follows the reading of the parameters file (now
502	from the \|stdin\|), the reading of the CT parameters file, and the
503	creation of the output file, where the processed data will be
504	stored.
505
506	*/
507
508	//++
509	// START
510	//--
511
512	// make unbuffered output
513
514	// cout.setf ( ios::stdio );
515
516	// parse command line options (see reflector.h)
517
518	parname[0] = '\0';
519
520	optarg = NULL;
521	while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
522	switch (ch) {
523	case 'f':
524	strcpy(parname, optarg);
525	break;
526	case 'h':
527	usage();
528	break;
529	default :
530	errflg++;
531	}
532
533	// show help if error
534
535	if ( errflg>0 )
536	usage();
537
538	// make some sort of presentation
539
540	present();
541
542	// read parameters file
543
544	if ( strlen(parname) < 1 )
545	readparam(NULL);
546	else
547	readparam(parname);
548
549	// read data from file or from STDIN?
550
551	Data_From_STDIN = get_data_from_stdin();
552
553	// get number of telescopes and camera geometries
554
555	ct_Number=get_ct_number();
556	if (ct_Number > MAX_NUMBER_OF_CTS)
557	{
558	error( SIGNATURE, "Number of telescopes is larger than maximum allowed (%i). Stoping camera program ...", MAX_NUMBER_OF_CTS);
559	exit(1);
560	}
561
562	for(int ict=0;ict<ct_Number;ict++){
563	ntrigger[ict]=0;
564	GeometryCamera[ict]=get_ct_geometry(ict);
565
566	//
567	// Get telescope coordinates (if supplied) from input card (units cm).
568	//
569	CTx[ict] = get_telescope_location_cm(ict,0);
570	CTy[ict] = get_telescope_location_cm(ict,1);
571	CTz[ict] = get_telescope_location_cm(ict,2);
572
573	//
574	// Get fraction of operative mirror:
575	//
576	mirror_frac[ict] = get_mirror_fraction(ict);
577	}
578
579	// structure holding the camera definition
580
581	TObjArray camgeom;
582
583	for (int i=0; i<ct_Number;i++){
584	switch(GeometryCamera[i]){
585	case 1: camgeom[i]=new MGeomCamMagic;
586	strcpy(GeometryName[i],"MGeomCamMagic");
587	TriggerPixels[i]=TRIGGER_PIXELS_1;
588	break;
589	case 2: camgeom[i]=new MGeomCamMagic919;
590	strcpy(GeometryName[i],"MGeomCamMagic919");
591	TriggerPixels[i]=TRIGGER_PIXELS_2;
592	break;
593	case 3: camgeom[i]=new MGeomCamMagicHG;
594	strcpy(GeometryName[i],"MGeomCamMagicHG");
595	TriggerPixels[i]=TRIGGER_PIXELS_3;
596	break;
597	case 5: camgeom[i]=new MGeomCamECO1000;
598	strcpy(GeometryName[i],"MGeomCamECO1000");
599	TriggerPixels[i]=TRIGGER_PIXELS_5;
600	break;
601	case 6: camgeom[i]=new MGeomCamECO1000HG;
602	strcpy(GeometryName[i],"MGeomCamECO1000HG");
603	TriggerPixels[i]=TRIGGER_PIXELS_6;
604	break;
605	case 8: camgeom[i]=new MGeomCamCT1;
606	strcpy(GeometryName[i],"MGeomCamCT1");
607	TriggerPixels[i]=TRIGGER_PIXELS_8;
608	break;
609	case 9: camgeom[i]=new MGeomCamCT1Daniel;
610	strcpy(GeometryName[i],"MGeomCamCT1Daniel");
611	TriggerPixels[i]=TRIGGER_PIXELS_9;
612	break;
613	default:
614	error( SIGNATURE, "Camera geometry %i is not defined. Stoping camera program ...", GeometryCamera[i]);
615	exit(1);
616	break;
617	}
618	}
619
620	ct_NPixels=0;
621
622	for(int ict=0;ict<ct_Number;ict++)
623	ct_NPixels=(((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels()>(UInt_t)ct_NPixels)?
624	(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels():
625	ct_NPixels;
626
627	// read WC and QE files
628
629	// FIX ME !
630	// Currently the PMT N of any camera with same average QE has the same QE.
631
632	QE = new float *** [ct_Number];
633	for(int ict=0;ict<ct_Number;ict++){
634	strcpy( qe_filename, get_qe_filename(ict));
635	read_QE(qe_filename,ict);
636	}
637
638	read_WC();
639
640	// write all images, even those without trigger?
641
642	Write_All_Event_Headers = get_write_all_event_headers();
643
644	Write_McEvt = get_write_McEvt() ;
645	Write_McTrig = get_write_McTrig() ;
646	Write_McFADC = get_write_McFadc() ;
647	Write_RawEvt = get_write_RawEvt() ;
648
649	FADC_Scan = get_FADC_Scan();
650	Trigger_Scan = get_Trigger_Scan();
651
652	Individual_Thres_Pixel = get_indi_thres_pixel();
653
654	get_secure_threhold(&riseDiskThres,&secureDiskThres);
655
656	get_FADC_properties
657	(&FADC_shape, &FADC_response_integ, &FADC_response_fwhm,
658	&FADC_shape_out, &FADC_resp_integ_out, &FADC_resp_fwhm_out,
659	&FADC_slices_per_ns, &FADC_slices_written);
660
661	// Allocate memory for the NSB contribution to the FADC signal:
662	nsb_fadcresp = new Float_t[(Int_t)(FADC_slices_per_ns*TOTAL_TRIGGER_TIME)];
663
664	FADC_high2low=get_High_to_Low();
665
666	// FIXME --- trigger properties may depend on the Camera geometry
667
668	get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);
669
670	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);
671
672	icontrigger =((int)((Trigger_loop_uthres-Trigger_loop_lthres)
673	/Trigger_loop_sthres)+1)*
674	(Trigger_loop_umult-Trigger_loop_lmult+1)*
675	(Trigger_loop_utop-Trigger_loop_ltop+1);
676
677	// Trigger loop operation is not implemented for Multi telescopes
678
679	if ( Trigger_Loop && ct_Number > 1 ){
680	cout<<"ERROR : camera::main : Trigger loop option is not";
681	cout<<" implemented for Multi Telescopes option. "<<icontrigger<<
682	" "<<ct_Number<<endl;
683	exit(1);
684	}
685
686	if (!Trigger_Loop){
687	get_Trigger_Single (qThreshold,
688	&Trigger_multiplicity[0],
689	&Trigger_topology[0]);
690	icontrigger=1;
691	}
692	else
693	get_threshold(qThreshold[0]);
694
695	// get filenames
696
697	if (! is_calibration_run())
698	for(int ict=0;ict<ct_Number;ict++)
699	{
700	strcpy( inname_CT[ict], get_input_filename(ict) );
701	if (strlen(inname_CT[ict]) == 0)
702	{
703	printf("\nError: missing input file name for CT id = %d. Exiting.\n\n", ict);
704	exit(1);
705	}
706	}
707
708	strcpy( starfieldname, get_starfield_filename() );
709	strcpy( datname, get_data_filename() );
710	strcpy( rootname, get_root_filename() );
711	strcpy( rootname_loop, get_loop_filename() );
712	strcpy( nsbpathname, get_nsb_directory() );
713	strcpy( nsbpath_outer, get_nsb_directory_outer() );
714
715	// get different parameters of the simulation
716
717	addElecNoise = add_elec_noise(&FADC_noise_inner, &FADC_noise_outer, &DIGITAL_noise, &Trigger_noise);
718	simulateNSB = get_nsb( &meanNSB, &nphe2NSB );
719	missPointing = get_misspointing(&missP_x,&missP_y);
720	Spotsigma = get_sigma_xy_cm_spot(&Spot_x,&Spot_y);
721	pmt_jitter = get_pmt_jitter_ns();
722
723
724	// get selections on the parameters
725
726	Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
727
728	//Parameters for starfield rotation
729
730	Starfield_rotate = get_starfield_rotate();
731
732	// log filenames information
733	for(Int_t ict=0;ict<ct_Number;ict++){
734
735	log(SIGNATURE,"\t%s : %i\n\t%20s:\t%s\n",
736	"------- TELESCOPE ",ict+1,
737	"Geometry ", GeometryName[ict]);
738	strcpy( qe_filename, get_qe_filename(ict));
739
740	printf("Telescope coordinates (cm): %.1f %.1f %.1f\n", CTx[ict], CTy[ict], CTz[ict]);
741
742	// Look to factor for NSB respect to emi_coat PMTs
743	if(strstr(qe_filename,"qe-emi-coat.RFL.dat") != 0)
744	factorqe_NSB[ict]=EMICOAT_NSB;
745	else if(strstr(qe_filename,"qe-emi.RFL.dat") != 0)
746	factorqe_NSB[ict]=EMI_NSB;
747	else if(strstr(qe_filename, "qe-intevac_hpd.RFL.dat") != 0)
748	factorqe_NSB[ict]=HPD_INTEVAC_NSB;
749	else if(strstr(qe_filename, "qe-hamamatsu_hpd.RFL.dat") != 0)
750	factorqe_NSB[ict]=HPD_HAMAMATSU_NSB;
751	else{
752	log(SIGNATURE,"\n\nWARNING!! : the factor for the diffuse NSB for this QE file (%s) is not known to the camera simulation. The number of phe(ns stated in the inputcard will be used.\n\n",qe_filename);
753	factorqe_NSB[ict]=1.0;
754	}
755	log(SIGNATURE,
756	"%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",
757	"Filenames",
758	"In", inname_CT[ict],
759	"Stars", starfieldname,
760	"NSB database","(inner pixels)", nsbpathname,
761	"(outer pixels)", nsbpath_outer,
762	"QE", qe_filename,
763	"Data", datname,
764	"ROOT", rootname
765	);
766
767	// log parameters information
768
769	log(SIGNATURE,
770	"%s:\n\t%20s: %f: %s\n\t%20s: %f\n\t%20s: %f\n",
771	"Parameters",
772	"NSB (phes/ ns 0.10.1 deg^2 239 m^2)", meanNSBfactorqe_NSB[ict], ONoff(simulateNSB),
773	"Pedestals = ", get_FADC_pedestal(),
774	"High to Low gain = ", FADC_high2low
775	);
776
777	// log Trigger information
778
779	if (Trigger_Loop) {
780	log(SIGNATURE,
781	"%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",
782	"Trigger Loop mode",
783	"Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres,
784	"Multiplicity",Trigger_loop_lmult,Trigger_loop_umult,
785	"Topology",Trigger_loop_ltop,Trigger_loop_utop,
786	rootname_loop);
787	}
788	else if (Individual_Thres_Pixel == FALSE){
789	log(SIGNATURE,
790	"%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n",
791	"Single Trigger mode",
792	"Threshold",qThreshold[ict][0],
793	"Multiplicity",Trigger_multiplicity[ict],
794	"Topology",Trigger_topology[ict]);
795	}
796	else{
797	log(SIGNATURE,
798	"%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n",
799	"Single Trigger mode",
800	"Threshold","Different for each pixel",
801	"Multiplicity",Trigger_multiplicity[ict],
802	"Topology",Trigger_topology[ict]);
803	}
804	log(SIGNATURE,"\t%s\n",
805	"END TELESCOPE --------");
806	}
807	// log flags information
808
809	log(SIGNATURE,
810	"%s:\n\t%23s: %s\n\t%23s: %s\n\t%23s: %s\n\t%23s: %s\n\t\t %s %3.2f,\n\t\t %s %3.2f (%s), %3.2f (%s) \n\t\t + %3.2f (%s)\n\n",
811	"Flags",
812	"Data_From_STDIN", ONoff(Data_From_STDIN),
813	"Write_All_Event_Headers", ONoff(Write_All_Event_Headers),
814	"Rotate Starfield", ONoff(Starfield_rotate),
815	"Electronic Noise", ONoff(addElecNoise),
816	"Trigger noise (amplitude relative to that of single phe): ",
817	Trigger_noise,
818	"FADC noise (ADC counts): ", FADC_noise_inner, "inner p.",
819	FADC_noise_outer, "outer p.",
820	DIGITAL_noise, "added digital noise at FADC"
821	);
822
823	if (! apply_gain_fluctuations())
824	{
825	log(SIGNATURE, "PMT gain fluctuations for the signal phes have been switched OFF");
826	log(SIGNATURE, "in the input card.\n\n");
827	}
828	if (! apply_noise_gain_fluctuations())
829	{
830	log(SIGNATURE, "PMT gain fluctuations for the NSB phes have been switched OFF\n");
831	log(SIGNATURE, "in the input card.\n\n");
832	}
833
834	if ( apply_gain_fluctuations() && !apply_noise_gain_fluctuations())
835	{
836	log(SIGNATURE, "Warning: you switched off PMT gain fluctuations for the NSB\n");
837	log(SIGNATURE, "photoelectrons, but not for the signal photoelectrons.\n\n");
838	}
839	else if ( !apply_gain_fluctuations() && apply_noise_gain_fluctuations())
840	{
841	log(SIGNATURE, "Warning: you switched off PMT gain fluctuations for the signal\n");
842	log(SIGNATURE, "photoelectrons, but not for the NSB photoelectrons.\n");
843	}
844
845
846	// log flags information
847
848	log(SIGNATURE,
849	"%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
850	"Root output",
851	"Write_McEvt", ONoff(Write_McEvt),
852	"Write_McTrig", ONoff(Write_McTrig),
853	"Write_McFADC", ONoff(Write_McFADC),
854	"Write_RawEvt", ONoff(Write_RawEvt));
855
856	// log selections
857
858	log(SIGNATURE,
859	"%s:\n\t%20s: %s (%f:%f)\n",
860	"Selections:",
861	"Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
862
863	// Definition and initialization of array to save trigger statistics
864
865	int ***ntriggerloop;
866
867	ntriggerloop= new int ** [(int)((Trigger_loop_uthres-Trigger_loop_lthres)
868	/Trigger_loop_sthres)+1];
869	for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;fthrescount+=Trigger_loop_sthres, ithrescount++){
870	ntriggerloop[ithrescount]= new int * [Trigger_loop_umult-Trigger_loop_lmult+1];
871	for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
872	ntriggerloop[ithrescount][imulticount]= new int [Trigger_loop_utop-Trigger_loop_ltop+1];
873	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
874	ntriggerloop[ithrescount][imulticount][itopocount]=0;
875	}
876	}
877	}
878
879	// We should be careful that topologies are sort from
880	// the less to the more restrictive one.
881
882	if (Trigger_loop_utop==Trigger_loop_ltop)
883	for(int is=0; is<3;is++)
884	isorttopo[is]=is;
885	else {
886	isorttopo[0]=1;
887	isorttopo[1]=0;
888	isorttopo[2]=2;
889	}
890
891	// get list of showers to evt. skip
892
893	nSkip = get_nskip_showers();
894
895	if (nSkip > 0) {
896	Skip = new int[ nSkip ];
897	get_skip_showers( Skip );
898
899	log(SIGNATURE, "There are some showers to skip:\n");
900	for (int i=0; i<nSkip; ++i)
901	log(SIGNATURE, "\tshower # %d\n", Skip[i]);
902	}
903
904
905	// read parameters from the ct.def file
906
907	//read_ct_file();
908
909	Int_t Lev0, Lev1;
910	Int_t Lev0MT[MAX_NUMBER_OF_CTS], Lev1MT[MAX_NUMBER_OF_CTS];
911
912	fadcValues = new TArrayC(FADC_slices_written);
913	fadcValuesLow = new TArrayC(FADC_slices_written);
914
915	// number of pixels for parameters
916
917	// Switched off writing TObject
918
919	MArray::Class()->IgnoreTObjectStreamer();
920	MParContainer::Class()->IgnoreTObjectStreamer();
921
922	// initialise ROOT
923
924	TROOT simple("simple", "MAGIC Telescope Monte Carlo");
925
926	// initialise instance of Trigger and FADC classes
927
928	MTrigger **Trigger_CT;
929	Trigger_CT = new MTrigger *[ct_Number];
930
931	for (int i=0; i<ct_Number;i++){
932	Trigger_CT[i] = new MTrigger(TriggerPixels[i],
933	((MGeomCam*)(camgeom.UncheckedAt(i))),
934	Trigger_gate_length,
935	Trigger_overlaping_time,
936	Trigger_response_ampl,
937	Trigger_response_fwhm, i); //@< A instance of the Class MTrigger
938	Trigger_CT[i]->SetSeed(UInt_t(i+get_seeds(0)));
939
940	if ( ! apply_gain_fluctuations())
941	Trigger_CT[i]->SetGainFluctuations(kFALSE);
942	}
943
944	for (int ict=0; ict<ct_Number;ict++){
945	// Generate database for trigger electronic noise
946	if (addElecNoise)
947	Trigger_CT[ict]->SetElecNoise(Trigger_noise) ;
948
949	// Set Right Discriminator threshold, taking into account trigger pixels
950	Trigger_CT[ict]->CheckThreshold(&qThreshold[ict][0],GeometryCamera[ict]);
951	}
952	// Set flag in pixel 0 (not used for trigger) that indicates if secure pixel
953	// is active: secureDiskThres*10000+riseDiskThres
954
955	for (int ict=0; ict<ct_Number;ict++){
956	if(riseDiskThres>0.0)
957	qThreshold[ict][0]=((UInt_t)secureDiskThres100)100+riseDiskThres;
958	else
959	qThreshold[ict][0]=0.0;
960	}
961	// Initialise McTrig information class if we want to save trigger informtion
962
963	MMcTrig **McTrig = NULL;
964	MMcTrigHeader **HeaderTrig = NULL;
965	MMcFadcHeader **HeaderFadc = NULL;
966
967	int numberBranches;
968
969	if (!Trigger_Loop)
970	numberBranches=ct_Number;
971	else
972	numberBranches=icontrigger;
973
974	if (Write_McTrig){
975
976	McTrig = new MMcTrig * [numberBranches];
977
978	for (int i=0;i<numberBranches;i++) {
979	McTrig[i] = new MMcTrig();
980	}
981
982	HeaderTrig = new MMcTrigHeader * [numberBranches];
983
984	for (int i=0;i<numberBranches;i++) {
985	HeaderTrig[i] = new MMcTrigHeader();
986	}
987	}
988
989	if (Write_McFADC){
990
991	HeaderFadc = new MMcFadcHeader * [numberBranches];
992
993	for (int i=0;i<numberBranches;i++) {
994	HeaderFadc[i] = new MMcFadcHeader();
995	}
996	}
997
998	float *fadc_jitter = new float[ct_Number];
999	for (int i=0; i<ct_Number;i++)
1000	fadc_jitter[i] = 0.;
1001
1002	MFadc **Fadc_CT;
1003	Fadc_CT = new MFadc *[ct_Number];
1004
1005	for (int i=0; i<ct_Number;i++)
1006	Fadc_CT[i] =
1007	new MFadc(((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels(),
1008	FADC_shape,
1009	FADC_response_integ,FADC_response_fwhm,
1010	FADC_shape_out,
1011	FADC_resp_integ_out,FADC_resp_fwhm_out,
1012	get_trig_delay(),
1013	FADC_slices_per_ns,
1014	FADC_slices_written); //@< A instance of the Class MFadc
1015
1016
1017	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1018	//
1019	// Set the FADC pedestals for that run
1020	// Some modifications
1021	// mut be done to simulate a more realistic distribution of the pedestals.
1022	// This simualtion is done int the SetPedestals methode inside the
1023	// class MFadc
1024	// Currentlly a input_pedestal array is declared with the pedestals.
1025	// Thy can also be set randomly following a flat distribution.
1026	//
1027	/////////////////////////////////////////////////////////////////////
1028
1029	Float_t input_pedestals[ct_NPixels];
1030
1031	for(int i=0;i<ct_NPixels;i++)
1032	input_pedestals[i]=get_FADC_pedestal();
1033	for (int ict=0; ict<ct_Number;ict++){
1034
1035	Fadc_CT[ict]->SetPedestals(input_pedestals);
1036	Fadc_CT[ict]->SetHigh2LowGain(FADC_high2low);
1037
1038	// Generate database for the Fadc electronic noise
1039
1040	if (!addElecNoise)
1041	continue;
1042
1043	MGeomCam camg = (MGeomCam)(camgeom.UncheckedAt(ict));
1044	UInt_t n_inner_pixels; // Number of inner(small) pixels
1045	for (n_inner_pixels = 0; n_inner_pixels < camg->GetNumPixels();
1046	n_inner_pixels ++)
1047	{
1048	if (camg->GetPixRatio(n_inner_pixels) < 1.)
1049	break;
1050	}
1051	Fadc_CT[ict]->SetElecNoise(FADC_noise_inner, FADC_noise_outer,
1052	n_inner_pixels);
1053	Fadc_CT[ict]->SetDigitalNoise(DIGITAL_noise);
1054	}
1055
1056	// Prepare the raw data output
1057
1058	// Header Tree
1059	MGeomCam *camdummy = new MGeomCam [ct_Number];
1060	for (int ict=0; ict<ct_Number;ict++){
1061	camdummy[ict]= ((MGeomCam*)(camgeom.UncheckedAt(ict)));
1062	}
1063	MRawRunHeader *RunHeader= new MRawRunHeader();
1064	MMcRunHeader *McRunHeader = new MMcRunHeader();
1065	MMcCorsikaRunHeader *McCorsikaRunHeader = new MMcCorsikaRunHeader("","",ct_Number);
1066
1067
1068	for (int i = 0; i < ct_Number; i++)
1069	McCorsikaRunHeader->FillCT(CTx[i], CTy[i], CTz[i], -1., -1., -1., -1., i);
1070
1071
1072	MMcConfigRunHeader **McConfigRunHeader = NULL;
1073	McConfigRunHeader = new MMcConfigRunHeader * [numberBranches];
1074	for (int i=0;i<numberBranches;i++) {
1075	McConfigRunHeader[i] = new MMcConfigRunHeader();
1076	}
1077
1078	// Header branch
1079
1080	MRawEvtHeader **EvtHeader = NULL;
1081
1082	if (Write_RawEvt) {
1083	EvtHeader = new MRawEvtHeader * [numberBranches];
1084
1085	for (int i=0;i<numberBranches;i++) {
1086	EvtHeader[i] = new MRawEvtHeader();
1087	}
1088
1089	}
1090
1091	// Data branch
1092
1093	MRawEvtData **EvtData = NULL; // Data branch
1094
1095	if (Write_RawEvt) {
1096	EvtData = new MRawEvtData * [numberBranches];
1097
1098	for (int i=0;i<numberBranches;i++) {
1099	EvtData[i] = new MRawEvtData();
1100	EvtData[i]->InitRead(RunHeader); // We need the RunHeader to read
1101	// number of pixels
1102	}
1103	}
1104
1105	MMcEvt **McEvt = NULL;
1106	MMcEvtBasic **McEvtBasic = NULL;
1107
1108	if (Write_McEvt)
1109	{
1110	McEvt = new MMcEvt *[ct_Number];
1111	McEvtBasic = new MMcEvtBasic *[ct_Number];
1112	for (int i = 0; i < ct_Number; i++)
1113	{
1114	McEvt[i] = new MMcEvt();
1115	McEvtBasic[i] = new MMcEvtBasic();
1116	}
1117	}
1118	//
1119	// initalize a temporal ROOT file
1120	//
1121
1122	TFile outfile_temp ( rootname , "RECREATE" );
1123
1124	// create a Tree for the Header Event
1125	TTree HeaderTree("RunHeaders","Headers of Run");
1126
1127	// define branches of Header Tree
1128
1129	char help[4];
1130
1131	HeaderTree.Branch("MRawRunHeader.","MRawRunHeader",
1132	&RunHeader);
1133
1134	HeaderTree.Branch("MMcRunHeader.","MMcRunHeader",
1135	&McRunHeader);
1136
1137	HeaderTree.Branch("MMcCorsikaRunHeader.","MMcCorsikaRunHeader",
1138	&McCorsikaRunHeader);
1139
1140	if(!Trigger_Loop && Write_McTrig && ct_Number==1){
1141
1142	HeaderTree.Branch("MMcTrigHeader.","MMcTrigHeader",
1143	&HeaderTrig[0]);
1144	}
1145
1146	if (ct_Number==1){
1147	// HeaderTree.Branch("MGeomCam.", "MGeomCam", &camdummy[0]);
1148	HeaderTree.Branch("MGeomCam.", GeometryName[0], &camdummy[0]);
1149	HeaderTree.Branch("MMcConfigRunHeader.","MMcConfigRunHeader",
1150	&McConfigRunHeader[0]);
1151	}
1152	else{
1153	char branchname[256];
1154	for (int ict=0; ict<ct_Number;ict++){
1155	sprintf(help,"%i",ict+1);
1156	strcpy (branchname, "MGeomCam;");
1157	strcat (branchname, & help[0]);
1158	strcat (branchname, ".");
1159	HeaderTree.Branch(branchname, GeometryName[ict],
1160	&camdummy[ict]);
1161	}
1162	for(int i=0;i<numberBranches;i++){
1163	sprintf(help,"%i",i+1);
1164	strcpy (branchname, "MMcConfigRunHeader;");
1165	strcat (branchname, & help[0]);
1166	strcat (branchname, ".");
1167	HeaderTree.Branch(branchname,"MMcConfigRunHeader",
1168	&McConfigRunHeader[i]);
1169	}
1170
1171	}
1172
1173	if ((Trigger_Loop \|\| ct_Number>1) && Write_McTrig){
1174	ibr=0;
1175	for(char branchname[256];ibr<numberBranches;ibr++){
1176
1177	sprintf(help,"%i",ibr+1);
1178	strcpy (branchname, "MMcTrigHeader;");
1179	strcat (branchname, & help[0]);
1180	strcat (branchname, ".");
1181	HeaderTree.Branch(branchname,"MMcTrigHeader",
1182	&HeaderTrig[ibr]);
1183	}
1184	}
1185
1186	if(!Trigger_Loop && Write_McFADC && ct_Number==1){
1187
1188	HeaderTree.Branch("MMcFadcHeader.","MMcFadcHeader",
1189	&HeaderFadc[0]);
1190	}
1191	if ((Trigger_Loop \|\| ct_Number>1) && Write_McFADC){
1192	ibr=0;
1193	for(char branchname[256];ibr<numberBranches;ibr++){
1194
1195	sprintf(help,"%i",ibr+1);
1196	strcpy (branchname, "MMcFadcHeader;");
1197	strcat (branchname, & help[0]);
1198	strcat (branchname, ".");
1199
1200	HeaderTree.Branch(branchname,"MMcFadcHeader",
1201	&HeaderFadc[ibr]);
1202	}
1203	}
1204
1205	// Fill branches for MRawRunHeader
1206
1207	RunHeader->SetMagicNumber(MRawRunHeader::kMagicNumber);
1208
1209	//
1210	// We set as format version 5, the current real data version as of February 2005.
1211	// These data have the peculiarity that the trigger pattern bits
1212	// (MRawEventHeader.fTrigPattern[0]) will be inverted (0 <-> 1) with respect
1213	// to the definitions in Mars/mtrigger/MTriggerPattern. We will account for this
1214	// below when setting the trigger pattern.
1215	//
1216	RunHeader->SetFormatVersion(5);
1217
1218
1219	RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
1220	RunHeader->SetRunType(256);
1221	RunHeader->SetRunNumber(0);
1222	RunHeader->SetNumSamples(FADC_slices_written, FADC_slices_written);
1223	RunHeader->SetNumCrates(1);
1224	RunHeader->SetNumPixInCrate(ct_NPixels);
1225
1226	// Fill branches for MMcTrigHeader
1227
1228	if(!Trigger_Loop && Write_McTrig && ct_Number==1){
1229
1230	HeaderTrig[0]->SetTopology((Short_t) Trigger_topology[0]);
1231	HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity[0]);
1232	HeaderTrig[0]->SetThreshold(qThreshold[0]);
1233	HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
1234	HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
1235	HeaderTrig[0]->SetOverlap(Trigger_overlaping_time);
1236	HeaderTrig[0]->SetGate(Trigger_gate_length);
1237	HeaderTrig[0]->SetElecNoise(Trigger_noise);
1238	HeaderTrig[0]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
1239	HeaderTrig[0]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
1240	}
1241
1242	if(!Trigger_Loop && Write_McTrig && ct_Number>1){
1243	for(int i=0;i<ct_Number;i++){
1244	HeaderTrig[i]->SetTopology((Short_t) Trigger_topology[i]);
1245	HeaderTrig[i]->SetMultiplicity((Short_t) Trigger_multiplicity[i]);
1246	HeaderTrig[i]->SetThreshold(qThreshold[i]);
1247	HeaderTrig[i]->SetAmplitud(Trigger_response_ampl);
1248	HeaderTrig[i]->SetFwhm(Trigger_response_fwhm);
1249	HeaderTrig[i]->SetOverlap(Trigger_overlaping_time);
1250	HeaderTrig[i]->SetGate(Trigger_gate_length);
1251	HeaderTrig[i]->SetElecNoise(Trigger_noise);
1252	HeaderTrig[i]->SetGainFluctuations(Trigger_CT[i]->GetGainFluctuations());
1253	HeaderTrig[i]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
1254	}
1255	}
1256	if(Trigger_Loop && Write_McTrig){
1257
1258	int iconcount;
1259	for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++,fthrescount+=Trigger_loop_sthres){
1260	for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
1261	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
1262	HeaderTrig[iconcount]->SetTopology((Short_t) isorttopo[itopocount+Trigger_loop_ltop]);
1263	HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
1264	for(int i=0;i<ct_NPixels;i++){
1265	fpixelthres[i]=
1266	((Float_t)(fthrescount)>=qThreshold[0][i])?
1267	(Float_t)(fthrescount):qThreshold[0][i];
1268	}
1269	HeaderTrig[iconcount]->SetThreshold( fpixelthres);
1270	HeaderTrig[iconcount]->SetAmplitud(Trigger_response_ampl);
1271	HeaderTrig[iconcount]->SetFwhm(Trigger_response_fwhm);
1272	HeaderTrig[iconcount]->SetOverlap(Trigger_overlaping_time);
1273	HeaderTrig[iconcount]->SetGate(Trigger_gate_length);
1274	HeaderTrig[iconcount]->SetElecNoise(Trigger_noise);
1275	HeaderTrig[iconcount]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
1276	HeaderTrig[iconcount]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
1277	iconcount++;
1278	}
1279	}
1280	}
1281	}
1282
1283	// Fill branches for MMcFadcHeader
1284	Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
1285
1286	if(!Trigger_Loop && Write_McFADC && ct_Number==1){
1287
1288	for(int k = 0; k < ct_NPixels; k++){
1289	if ( ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetPixRatio(k) < 1.)
1290	fadc_elecnoise[k]=FADC_noise_outer; // outer pixels
1291	else
1292	fadc_elecnoise[k]=FADC_noise_inner; // inner pixels
1293
1294	fadc_diginoise[k]=DIGITAL_noise;
1295	}
1296
1297	HeaderFadc[0]->SetShape((Float_t)FADC_shape);
1298	HeaderFadc[0]->SetShapeOuter((Float_t)FADC_shape_out);
1299	HeaderFadc[0]->SetAmplitud(FADC_response_integ,
1300	FADC_resp_integ_out);
1301	HeaderFadc[0]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
1302	HeaderFadc[0]->SetLow2High(FADC_high2low);
1303	HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],
1304	((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
1305	HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
1306	((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
1307
1308	//
1309	// Fill also the flag indicating whether the PMT gain fluctuations
1310	// were simulated or not:
1311	//
1312	HeaderFadc[0]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
1313	HeaderFadc[0]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
1314
1315	}
1316
1317	if(!Trigger_Loop && Write_McFADC && ct_Number>1){
1318	for(int i=0;i<ct_Number;i++){
1319
1320	for(int k = 0; k < ct_NPixels; k++){
1321	if ( ((MGeomCam*)(camgeom.UncheckedAt(i)))->GetPixRatio(k) < 1.)
1322	fadc_elecnoise[k]=FADC_noise_outer; // outer pixels
1323	else
1324	fadc_elecnoise[k]=FADC_noise_inner; // inner pixels
1325
1326	fadc_diginoise[k]=DIGITAL_noise;
1327	}
1328
1329	Fadc_CT[i]->GetPedestals(&fadc_pedestals[0]);
1330	HeaderFadc[i]->SetShape((Float_t)FADC_shape);
1331	HeaderFadc[i]->SetShapeOuter((Float_t)FADC_shape_out);
1332	HeaderFadc[i]->SetAmplitud(FADC_response_integ,
1333	FADC_resp_integ_out);
1334	HeaderFadc[i]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
1335	HeaderFadc[i]->SetLow2High(FADC_high2low);
1336	HeaderFadc[i]->SetPedestal(&fadc_pedestals[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
1337	HeaderFadc[i]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
1338	((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
1339
1340	//
1341	// Fill also the flag indicating whether the PMT gain fluctuations
1342	// were simulated or not:
1343	//
1344	HeaderFadc[i]->SetGainFluctuations(Trigger_CT[i]->GetGainFluctuations());
1345	HeaderFadc[i]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
1346
1347	}
1348	}
1349
1350	if(Trigger_Loop && Write_McFADC){
1351
1352	for(int k = 0; k < ct_NPixels; k++){
1353	if ( ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetPixRatio(k) < 1.)
1354	fadc_elecnoise[k]=FADC_noise_outer; // outer
1355	else
1356	fadc_elecnoise[k]=FADC_noise_inner; // inner
1357
1358	fadc_diginoise[k]=DIGITAL_noise;
1359	}
1360
1361	int iconcount;
1362	for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1363	for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
1364	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
1365	Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
1366	HeaderFadc[iconcount]->SetShape((Float_t)FADC_shape);
1367	HeaderFadc[iconcount]->SetShapeOuter((Float_t)FADC_shape_out);
1368	HeaderFadc[iconcount]->SetAmplitud(FADC_response_integ,
1369	FADC_resp_integ_out);
1370	HeaderFadc[iconcount]->SetFwhm(FADC_response_fwhm,
1371	FADC_resp_fwhm_out);
1372	HeaderFadc[iconcount]->SetLow2High(FADC_high2low);
1373	HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0], ct_NPixels);
1374	HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],
1375	&fadc_diginoise[0],ct_NPixels);
1376
1377	HeaderFadc[iconcount]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
1378	HeaderFadc[iconcount]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
1379
1380	iconcount++;
1381	}
1382	}
1383	}
1384	}
1385
1386	// UnSet flag in pixel 0 (not used for trigger) that indicates if secure pixel
1387	// is active once it has been stored
1388
1389	for (int ict=0; ict<ct_Number;ict++){
1390	qThreshold[ict][0]=999999.99;
1391	}
1392
1393	// create a Tree for the Event data stream
1394	TTree EvtTree("Events","Normal Triggered Events");
1395
1396	if (Write_McEvt && ct_Number == 1)
1397	{
1398	EvtTree.Branch("MMcEvt.","MMcEvt", &McEvt[0]);
1399	EvtTree.Branch("MMcEvtBasic.","MMcEvtBasic", &McEvtBasic[0]);
1400	}
1401
1402	if (Write_McEvt && ct_Number!=1)
1403	{
1404	ibr=0;
1405	for(char branchname[256]; ibr < numberBranches; ibr++)
1406	{
1407	sprintf(help,"%i",ibr+1);
1408	strcpy (branchname, "MMcEvt;");
1409	strcat (branchname, & help[0]);
1410	strcat (branchname, ".");
1411	EvtTree.Branch(branchname,"MMcEvt",
1412	&McEvt[ibr]);
1413	}
1414	for(char branchname[256]; ibr < numberBranches; ibr++)
1415	{
1416	sprintf(help,"%i",ibr+1);
1417	strcpy (branchname, "MMcEvtBasic;");
1418	strcat (branchname, & help[0]);
1419	strcat (branchname, ".");
1420	EvtTree.Branch(branchname,"MMcEvtBasic",
1421	&McEvtBasic[ibr]);
1422	}
1423	}
1424
1425	if(!Trigger_Loop && ct_Number==1){
1426
1427	if (Write_RawEvt){
1428	EvtTree.Branch("MRawEvtHeader.","MRawEvtHeader",
1429	&EvtHeader[0]);
1430	EvtTree.Branch("MRawEvtData.","MRawEvtData",
1431	&EvtData[0]);
1432	}
1433	if (Write_McTrig){
1434	EvtTree.Branch("MMcTrig.","MMcTrig",
1435	&McTrig[0]);
1436	}
1437	}
1438	else{
1439
1440	if (Write_McTrig){
1441	ibr=0;
1442	for(char branchname[256];ibr<numberBranches;ibr++){
1443
1444	sprintf(help,"%i",ibr+1);
1445	strcpy (branchname, "MMcTrig;");
1446	strcat (branchname, & help[0]);
1447	strcat (branchname, ".");
1448	EvtTree.Branch(branchname,"MMcTrig",
1449	&McTrig[ibr]);
1450	}
1451	}
1452	}
1453
1454	if ((Trigger_Loop \|\| ct_Number>1) && Write_RawEvt){
1455	ibr=0;
1456	for(char branchname[256];ibr<numberBranches;ibr++){
1457
1458	sprintf(help,"%i",ibr+1);
1459	strcpy (branchname, "MRawEvtHeader;");
1460	strcat (branchname, & help[0]);
1461	strcat (branchname, ".");
1462	EvtTree.Branch(branchname,"MRawEvtHeader",
1463	&EvtHeader[ibr]);
1464	}
1465	ibr=0;
1466	for(char branchname[256];ibr<numberBranches;ibr++){
1467
1468	sprintf(help,"%i",ibr+1);
1469	strcpy (branchname, "MRawEvtData;");
1470	strcat (branchname, & help[0]);
1471	strcat (branchname, ".");
1472	EvtTree.Branch(branchname,"MRawEvtData",
1473	&EvtData[ibr]);
1474	}
1475	}
1476
1477
1478	TApplication theAppTrigger("App", &argc, argv);
1479
1480	if(FADC_Scan){
1481	if (gROOT->IsBatch()) {
1482	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1483	// return 1;
1484	}
1485	}
1486	if(FADC_Scan){
1487	//TApplication theAppFadc("App", &argc, argv);
1488
1489	if (gROOT->IsBatch()) {
1490	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1491	// return 1;
1492	}
1493	}
1494
1495	// prepare the NSB simulation
1496
1497	// Instance of the Mlons class
1498
1499	MLons lons(0, Trigger_response_ampl, Trigger_response_fwhm,
1500	FADC_shape, FADC_response_integ, FADC_response_fwhm,
1501	FADC_slices_per_ns, apply_noise_gain_fluctuations());
1502
1503	lons.SetSeed(((UInt_t)(get_seeds(1)get_seeds(1)get_seeds(0))));
1504
1505	lons.SetPath(nsbpathname);
1506
1507	// Instance of the Mlons class
1508	MLons lons_outer(0, Trigger_response_ampl, Trigger_response_fwhm,
1509	FADC_shape_out,FADC_resp_integ_out,FADC_resp_fwhm_out,
1510	FADC_slices_per_ns, apply_noise_gain_fluctuations());
1511
1512	lons_outer.SetSeed(((UInt_t)(get_seeds(1)get_seeds(0)get_seeds(0))));
1513
1514	lons_outer.SetPath(nsbpath_outer);
1515
1516	if( simulateNSB){
1517
1518	//
1519	// Calculate the non-diffuse NSB photoelectron rates
1520	//
1521
1522	//
1523	// FIXME! --- star NSB different for each camera?
1524	// Then we will have to use mirror_frac[ict]
1525	//
1526	log(SIGNATURE,"Produce NSB rates from Star Field...\n");
1527
1528	k = produce_nsbrates( starfieldname,
1529	((MGeomCam*)(camgeom.UncheckedAt(0))),
1530	nsbrate_phepns,
1531	0,
1532	mirror_frac[0]);
1533
1534	//
1535	// Call to "produce_nsbrates" above accounts ONLY for
1536	// non-diffuse NSB (i.e. from stars). NOTE: produce_nsbrates already
1537	// accounts for the possibly different light collection efficiencies
1538	// of inner and outer pixels, through a call (see function) to
1539	// "produce_phes". The output array nsbrate_phepns contains only the
1540	// rates due to individual stars in the FOV, and not diffuse NSB light!
1541	//
1542
1543	if (k != 0){
1544	cout << "Error when reading starfield... \nExiting.\n";
1545	exit(1);
1546	}
1547
1548	// calculate diffuse rate correcting for the pixel size and telescope
1549
1550	float factorNSB_ct;
1551	for(int ict=0;ict<ct_Number;ict++){
1552
1553	// First we set the factor due to the mirror size with respect to the normal
1554	// MAGIC mirror geometry.
1555
1556	switch(GeometryCamera[ict]){
1557	case 1:
1558	case 2:
1559	case 3:
1560	case 4:
1561	factorNSB_ct=1.0;
1562	break;
1563	case 5:
1564	case 6:
1565	case 7:
1566	factorNSB_ct=1000.0/239.0;
1567	break;
1568	case 8:
1569	case 9:
1570	default:
1571	factorNSB_ct=1.0;
1572	break;
1573	}
1574
1575	factorNSB_ct = factorNSB_ct /
1576	((((MGeomCam)(camgeom.UncheckedAt(ict)))).GetCameraDist()1000
1577	(((MGeomCam)(camgeom.UncheckedAt(ict)))).GetCameraDist()1000
1578	PIXEL_SIZE*PIXEL_SIZE); // [mm^-2]
1579
1580	for(UInt_t ui=0;
1581	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); ui++){
1582	const Float_t size=
1583	(((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD();
1584	// Returns distance [mm] between two paralel sides of pixel
1585
1586
1587	Float_t diffusensb = meanNSB*mirror_frac[ict];
1588
1589	//
1590	// If pixel is outer pixel:
1591	//
1592	if ( ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetPixRatio(ui) < 1. )
1593	diffusensb *= WC_outer[1][90] / WC[1][90];
1594	//
1595	// FIXME! Correction above is for (possibly) different light collection efficiencies of
1596	// inner and outer pixels. For the moment we assume the angular dependence of the light
1597	// collection is the same and hence use simply the ratio of efficiencies for light impinging
1598	// perpendicular to the camera plane (index 90 stands for 90 degrees)
1599	//
1600
1601	diffnsb_phepns[ict][ui] = (Int_t(diffusensbfactorNSB_ct100sizesize+0.5))/(100.0)*factorqe_NSB[ict];
1602	}
1603	}
1604
1605	// calculate nsb rate including diffuse and starlight
1606	// we also include the extinction effect
1607	for(int ict=0;ict<ct_Number;ict++){
1608	for(j=0;j<iNUMWAVEBANDS;j++){
1609	// calculate the effect of the atmospheric extinction (for stars!)
1610
1611	zenfactor = pow(10., -0.4 * ext[j] );
1612
1613	for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++)
1614	{
1615	nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS +
1616	zenfactor * nsbrate_phepns[ui][j];
1617	nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
1618	}
1619	}
1620	}
1621	}
1622
1623
1624	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1625	//
1626	// Now 500 empty events with the condition in which the camera is run
1627	// are simulated. In this way one gets an estimation of the
1628	// fluctuations of the pedestal of each FADC channel.
1629	// This computation is done assuming any noise that affects
1630	// the FADC but there is no rotation of the Star Field (otherwise it
1631	// should be done for each event). So it is valid if no starfield
1632	// rotation is used.
1633	//
1634	// Changed 20/03/2004, AM: now we no longer calculate the individual
1635	// FADC slice RMS. Due to correlations in the noise of neighboring
1636	// slides, it follows that RMS(sum_n_slices) != sqrt(n)*RMS(single_slice)
1637	//
1638	// In the analysis in Mars, however, the RMS of the fluctuations of the
1639	// signal (resulting from the integration of n slices) is estimated as
1640	// sqrt(n) * MMcFadcHeader.fPedesSigmaHigh, where the latter value,
1641	// stored in the camera output, is calculated in the next lines of code.
1642	// We have then made the following, as is being done also in real data:
1643	// calculate the RMS of the distribution of the sum of 14 slices, then
1644	// the stored value is
1645	// MMcFadcHeader.fPedesSigmaHigh = RMS(sum_14_slices)/sqrt(14), and the
1646	// same for the low gain. It can be seen that the RMS of the sum of n
1647	// slices, with n not too low (n>=6 or so), is more or less:
1648	// RMS(sum_n_slices) ~ sqrt(n) * RMS(sum_14_slices)/sqrt(14)
1649	//
1650	// The reason to sum 14 slices (and not 15) comes from the fact that in
1651	// real data there is a FADC clock noise affecting differently odd and
1652	// even FADC slices, so always an even number of them is added up so that
1653	// this cancels out. So we do the calculation in the same way as in real
1654	// data.
1655	//
1656	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1657
1658	Int_t empty_events = 500;
1659	Int_t n_slices = 14;
1660
1661	for(int ict=0;ict<ct_Number;ict++){
1662	for(UInt_t ui=0;
1663	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1664	ui++){
1665	fadc_sigma[ui]=0.0;
1666	fadc_sigma_low[ui]=0.0;
1667	}
1668	for(int ie=0; ie < empty_events; ie++){
1669	Fadc_CT[ict]->Reset();
1670	if (addElecNoise){
1671	Fadc_CT[ict]->ElecNoise() ;
1672	}
1673	if(simulateNSB)
1674	{
1675	for(UInt_t ui=0;
1676	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1677	ui++)
1678	{
1679	if(nsb_phepns[ict][ui]>0.0)
1680	{
1681	if((((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD() >
1682	(((MGeomCam)(camgeom.UncheckedAt(ict))))[0].GetD())
1683	{
1684
1685	k = lons_outer.GetResponse(nsb_phepns[ict][ui],0.01,
1686	& nsb_trigresp[0],
1687	& nsb_fadcresp[0]);
1688	}
1689	else
1690	{
1691	k = lons.GetResponse(nsb_phepns[ict][ui],0.01,
1692	& nsb_trigresp[0],& nsb_fadcresp[0]);
1693	}
1694
1695	if(k==0)
1696	{
1697	cout << "Exiting.\n";
1698	exit(1);
1699	}
1700	Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
1701	}
1702	}
1703	}
1704	Fadc_CT[ict]->Pedestals();
1705
1706	for(UInt_t ui=0;
1707	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1708	ui++)
1709	{
1710	//
1711	// We add up n_slices FADC slices (pedestal subtracted), then
1712	// calculate the sigma_n-1 of this sum for the number of generated
1713	// noise events (=empty_events).
1714	//
1715	Float_t sumslices = Fadc_CT[ict]->AddNoiseInSlices(ui,1,n_slices);
1716	fadc_sigma[ui] += sumslices*sumslices;
1717
1718	// Now the low gain:
1719	sumslices = Fadc_CT[ict]->AddNoiseInSlices(ui,0,n_slices);
1720	fadc_sigma_low[ui]+= sumslices*sumslices;
1721	}
1722	}
1723
1724	for(UInt_t ui=0;
1725	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1726	ui++)
1727	{
1728	Float_t s_high = fadc_sigma[ui] / (Float_t)(empty_events-1) /
1729	(Float_t)(n_slices);
1730	fadc_sigma[ui] = sqrt(s_high);
1731
1732	Float_t s_low = fadc_sigma_low[ui] / (Float_t)(empty_events-1) /
1733	(Float_t)(n_slices);
1734	fadc_sigma_low[ui] = sqrt(s_low);
1735	}
1736
1737	HeaderFadc[ict]->SetPedestalSigma(&fadc_sigma_low[0],&fadc_sigma[0],
1738	((MGeomCam*)(camgeom.UncheckedAt(ict)))
1739	->GetNumPixels());
1740	}
1741
1742
1743	if (is_calibration_run())
1744	{
1745	DoCalibration(Fadc_CT, Trigger_CT, camgeom, nsb_trigresp, nsb_fadcresp,
1746	&lons, &lons_outer, nsb_phepns, addElecNoise,
1747	&EvtTree, EvtHeader, McEvt, EvtData);
1748
1749	HeaderTree.Fill() ;
1750
1751	outfile_temp.Write();
1752	outfile_temp.Close();
1753
1754	cout << endl << "Calibration run finished!" << endl << endl;
1755
1756	return (0);
1757	}
1758
1759	//
1760	// Read the reflector file with the Cherenkov data
1761	//
1762
1763	// select input file
1764
1765	if ( Data_From_STDIN )
1766	inputfile[0] = stdin;
1767
1768	else
1769	{
1770	for(int ict = 0; ict < ct_Number; ict++)
1771	{
1772	log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
1773	inputfile[ict] = fopen( inname_CT[ict], "r" );
1774
1775	if ( inputfile[ict] == NULL )
1776	error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
1777	}
1778	}
1779
1780	// get signature, and check it
1781
1782	for(int ict=0;ict<ct_Number;ict++)
1783	{
1784	if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE)
1785	exit(1);
1786	}
1787
1788	// open data file
1789
1790	log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
1791	datafile.open( datname );
1792
1793	if ( datafile.bad() )
1794	error( SIGNATURE, "Cannot open data file: %s\n", datname );
1795
1796	// initializes flag
1797
1798	strcpy( flag, " \0" );
1799	strcpy( flag_new, " \0" );
1800
1801	// allocate space for PMTs numbers of pixels
1802
1803	fnpix = new float [ct_NPixels];
1804
1805
1806	//!@' @#### Main loop.
1807	//@'
1808
1809	// get flag
1810
1811	for(int ict=0;ict<ct_Number;ict++)
1812	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1813
1814
1815	// loop over the file
1816
1817	still_in_loop = TRUE;
1818
1819	// FIXME --- check if not eof for all input reflector files
1820
1821	while (
1822	((! Data_From_STDIN) && ( !feof(inputfile[0]) ))
1823	\|\|
1824	(Data_From_STDIN && still_in_loop)
1825	) {
1826
1827	// reading .rfl files
1828	if(!isA( flag, FLAG_START_OF_RUN )){
1829
1830	// We break the main loop
1831	cout<<"Warning: Expected start of run flag, but found:"<<flag<<endl;
1832	cout<<" We break the main loop"<<endl;
1833	break;
1834	}
1835	else { // found start of run
1836
1837	for(int ict=0;ict<ct_Number;ict++) {
1838
1839	fread( flag_new, 4, 1, inputfile[ict] );
1840
1841	if(!isA( flag_new, FLAG_START_OF_HEADER)){
1842
1843	// We break the main loop
1844	cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
1845	cout<<" We break the main loop"<<endl;
1846	break;
1847	}
1848
1849	// fread((char)&mcrunh,(SIZE_OF_MCRUNHEADER)sizeof(float), 1, inputfile[ict] );
1850	// AM, changed reading method, 02/2004:
1851
1852	mcrunh.read(inputfile[ict]);
1853
1854	}
1855
1856	nshow=0;
1857
1858	for(int ict=0;ict<ct_Number;ict++){
1859	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1860	}
1861	while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
1862	for(int ict=0;ict<ct_Number;ict++) {
1863	fread( flag_new, 4, 1, inputfile[ict] );
1864	if(!isA( flag_new, FLAG_EVENT_HEADER)){
1865
1866	// We break while events loop
1867	cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<" "<<ict<<endl;
1868	cout<<" We break the while events loop"<<endl;
1869	break;
1870	}
1871	}
1872
1873
1874	//
1875	// Clear Trigger and Fadc
1876	//
1877	for(int ict=0;ict<ct_Number;ict++) {
1878	Trigger_CT[ict]->Reset() ;
1879	Trigger_CT[ict]->ClearFirst();
1880	Trigger_CT[ict]->ClearZero();
1881	Fadc_CT[ict]->Reset() ;
1882	}
1883
1884	++nshow;
1885	if((nshow+ntshow+1)%100 == 1)
1886	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
1887
1888	// get MCEventHeader
1889	if (reflector_file_version<6){
1890	for(int ict=0;ict<ct_Number;ict++)
1891	{
1892	mcevth[ict].read(inputfile[ict]);
1893	}
1894
1895	}
1896	else{
1897	for(int ict=0;ict<ct_Number;ict++)
1898	{
1899	mcevth_2[ict].read(inputfile[ict]);
1900	}
1901	}
1902
1903	//
1904	// AM March 2004 simplified impact parameter calculation,
1905	// and allowed for correct estimate also for telescopes not
1906	// placed at Corsika's origin of coordinates (0.,0.).
1907	// FIXME: telescope coordinates still set by the user in the
1908	// input card, since they are not available in the reflector
1909	// file!
1910	//
1911	// Calculate impact parameter as distance between telescope
1912	// location and shower axis. In the previous implementation,
1913	// the definition was the distance between telescope axis and
1914	// shower axis, but this has less physical meaning in general.
1915	// Light yield depends above all on distance to shower axis!
1916	// Of course for shower axis paralel to telescope the old and
1917	// the new definitions of impact parameter agree.
1918	//
1919
1920
1921	// read the direction of the incoming shower
1922	// It is done only for one telescope since it is suposed
1923	// to be the same shower for all of them
1924	if (reflector_file_version<6){
1925	thetashw = mcevth[0].get_theta();
1926	phishw = mcevth[0].get_phi();
1927	}
1928	else{
1929	thetashw = mcevth_2[0].get_theta();
1930	phishw = mcevth_2[0].get_phi();
1931	}
1932	Zenith=thetashw; Azimutal=phishw;
1933
1934	// calculate vector for shower
1935
1936	l1 = sin(thetashw)*cos(phishw);
1937	m1 = sin(thetashw)*sin(phishw);
1938	n1 = cos(thetashw);
1939
1940	//
1941	// Note, A.M.:
1942	// Attention! "core vector from mcevth*.get_core
1943	// is the vector from the core to the telescope!
1944	//
1945	Float_t core2ct_x;
1946	Float_t core2ct_y;
1947
1948	if (reflector_file_version<6)
1949	mcevth[0].get_core(&core2ct_x, &core2ct_y);
1950	else
1951	mcevth_2[0].get_core(&core2ct_x, &core2ct_y);
1952
1953	//
1954	// Then true core position in Corsika's system is:
1955	//
1956
1957	coreX = CTx[0] - core2ct_x;
1958	coreY = CTy[0] - core2ct_y;
1959
1960	//
1961	// FIXME: This may not work fine for CTs at z != 0 !!
1962	//
1963	for(int ict=0;ict<ct_Number;ict++)
1964	impactD[ict] = dist_r_P( CTx[ict], CTy[ict], CTz[ict],
1965	l1, m1, n1, coreX, coreY, 0. );
1966
1967
1968	// energy cut
1969
1970	if ( Select_Energy ) {
1971	if (reflector_file_version<6)
1972	if (( mcevth[0].get_energy() < Select_Energy_le ) \|\|
1973	( mcevth[0].get_energy() > Select_Energy_ue )) {
1974	log(SIGNATURE, "select_energy: shower rejected.\n");
1975	continue;
1976	}
1977	else
1978	if (( mcevth_2[0].get_energy() < Select_Energy_le ) \|\|
1979	( mcevth_2[0].get_energy() > Select_Energy_ue )) {
1980	log(SIGNATURE, "select_energy: shower rejected.\n");
1981	continue;
1982	}
1983	}
1984
1985	inumphe=0;
1986
1987	for(int ict=0;ict<ct_Number;ict++){
1988
1989	// Read first and last time and put inumphe_CT[0] to 0
1990
1991	if (reflector_file_version<6)
1992	mcevth[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1993	else
1994	mcevth_2[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1995
1996	inumphe_CT[ict]=0;
1997
1998
1999	// Obtain the FADC jitter of 1 FADC slice. This is a time to be added to the
2000	// time of all photons in an event, before digitalization of the signal. It is
2001	// therefore the same time shift for all pixels in a CT.
2002
2003	fadc_jitter[ict] =
2004	(1./Fadc_CT[ict]->GetFadcSlicesPerNanosec()) * RandomNumber; //ns
2005
2006	// read the photons and produce the photoelectrons
2007
2008	k = produce_phes( inputfile[ict],
2009	((MGeomCam*)(camgeom.UncheckedAt(ict))),
2010	WAVEBANDBOUND1,
2011	WAVEBANDBOUND6,
2012	Trigger_CT[ict], // will be changed by the function!
2013	Fadc_CT[ict], // will be changed by the function!
2014	&inumphe_CT[ict], // important for later: the size of photoe[]
2015	fnpix, // will be changed by the function!
2016	&ncph[ict], // will be changed by the function!
2017	&arrtmin_ns, // will be changed by the function!
2018	&arrtmax_ns, // will be changed by the function!
2019	ict,
2020	mirror_frac[ict],
2021	fadc_jitter[ict]);
2022
2023	inumphe = (inumphe < inumphe_CT[ict])? inumphe_CT[ict] : inumphe;
2024	// inumphe will be the maximum number of phe in any of the telescopes
2025
2026
2027	if( k != 0 ){ // non-zero return value means error
2028	cout << "Exiting.\n";
2029	exit(1);
2030	}
2031	}
2032
2033	// NSB simulation
2034
2035	if(simulateNSB && inumphe >= nphe2NSB)
2036	{
2037
2038	if(Starfield_rotate){
2039
2040	// Introduction rho angle
2041
2042	zenith = thetashw;
2043	azimutal = phishw;
2044	C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
2045	C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
2046	C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) +
2047	C3 * C3 );
2048	rho = acos( C1/C2 );
2049
2050	if ( sin(azimutal) < 0)
2051	rho = 2 * 3.14159 - rho;
2052	else
2053	rho = rho;
2054
2055	rho = rho*180/3.14159;
2056
2057	// Rotation of the NSB
2058	// FIXME --- We should rotate for all cameras. Is it always the same rho?
2059	for(int ict=0;ict<ct_Number;ict++)
2060	k = size_rotated(&nsb_phepns_rotated[ict][0],
2061	nsb_phepns[ict],
2062	rho);
2063	}
2064
2065	// Fill trigger and fadc response in the trigger class from the database
2066	for(int ict=0;ict<ct_Number;ict++)
2067	{
2068	for(UInt_t ui=0;
2069	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
2070	ui++)
2071	{
2072	if(nsb_phepns_rotated[ict][ui]>0.0)
2073	{
2074	if((((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD() >
2075	(((MGeomCam)(camgeom.UncheckedAt(ict))))[0].GetD())
2076	{
2077	k=lons_outer.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
2078	& nsb_trigresp[0],
2079	& nsb_fadcresp[0]);
2080	}
2081	else
2082	{
2083	k=lons.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
2084	& nsb_trigresp[0],& nsb_fadcresp[0]);
2085	}
2086	if(k==0)
2087	{
2088	cout << "Exiting.\n";
2089	exit(1);
2090	}
2091	Trigger_CT[ict]->AddNSB(ui,nsb_trigresp);
2092	Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
2093	}
2094	}
2095	}
2096
2097	}// end if(simulateNSB && inumphe >= nphe2NSB) ...
2098
2099
2100	for(int ict=0;ict<ct_Number;ict++)
2101	{
2102	inumphensb[ict]=0;
2103
2104	for (UInt_t ui=0;
2105	ui < ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
2106	ui++)
2107	inumphensb[ict]+=nsb_phepns[ict][ui]*TOTAL_TRIGGER_TIME;
2108
2109	ntcph[ict]+=ncph[ict];
2110	if ((nshow+ntshow+1)%100 == 1){
2111	log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
2112	ncph[ict], ntcph[ict]);
2113
2114	cout << "Total number of phes in CT "<<ict<<": "
2115	<< inumphe_CT[ict]<<" (+ ";
2116	cout<<inumphensb[ict]<<" mean expected number from NSB)"<<endl;
2117	}
2118	}
2119
2120	// skip it ?
2121
2122	int i;
2123	for (i=0; i<nSkip; ++i ) {
2124	if (Skip[i] == (nshow+ntshow)) {
2125	i = -1;
2126	break;
2127	}
2128	}
2129
2130	// if after the previous loop, the exit value of i is -1
2131	// then the shower number is in the list of showers to be
2132	// skipped
2133
2134	if (i == -1) {
2135	log(SIGNATURE, "\t\tskipped!\n");
2136	continue;
2137	}
2138
2139	//++++++++++++++++++++++++++++++++++++++++++++++++++
2140	// at this point we have a camera full of
2141	// ph.e.s
2142	// we should first apply the trigger condition,
2143	// and if there's trigger, then clean the image,
2144	// calculate the islands statistics and the
2145	// other parameters of the image (Hillas' parameters
2146	// and so on).
2147	//--------------------------------------------------
2148
2149	// TRIGGER HERE
2150
2151	// We should simulate the AC coupling behaviour:
2152	// For the FADC it is only done for the NSB while producing
2153	// the StarResponse database.
2154	// For the trigger is done in the Trigger.Diskriminate(), which
2155	// is called later (it should be separated to speed up the program)
2156	//
2157
2158	// now the noise of the electronic
2159	// (preamps, optical transmission,..) is introduced.
2160	// This is done inside the class MTrigger by the method ElecNoise.
2161	//
2162
2163	for(int ict=0;ict<ct_Number;ict++)
2164	{
2165	if (addElecNoise && inumphe >= nphe2NSB)
2166	{
2167	Trigger_CT[ict]->ElecNoise(Trigger_noise);
2168	Fadc_CT[ict]->ElecNoise();
2169	}
2170	}
2171
2172	// now a shift in the fadc signal due to the pedestals is
2173	// introduced
2174	// This is done inside the class MFadc by the method Pedestals
2175
2176	for(int ict=0;ict<ct_Number;ict++)
2177	Fadc_CT[ict]->Pedestals();
2178
2179
2180	// We study several trigger conditons
2181	if(Trigger_Loop)
2182	{
2183	// Set to zero the flag to know if some conditon has triggered
2184	btrigger=0;
2185	flagstoring = 0;
2186
2187	// Loop over trigger threshold
2188	int iconcount;
2189	for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;
2190	fthrescount <= Trigger_loop_uthres;
2191	ithrescount++, fthrescount += Trigger_loop_sthres)
2192	{
2193	for (int i=0;i<ct_NPixels;i++)
2194	{
2195	fpixelthres[i] =
2196	((Float_t)(fthrescount)>=qThreshold[0][i])?
2197	(Float_t)(fthrescount):qThreshold[0][i];
2198
2199	// Rise the discrimnator threshold to avoid huge rates
2200
2201	if(riseDiskThres>0.0 && simulateNSB && inumphe >= nphe2NSB)
2202	for(int ii=0;ii<ct_NPixels;ii++)
2203	if( nsb_phepns_rotated[0][ii]>riseDiskThres)
2204	fpixelthres[ii]=secureDiskThres;
2205	}
2206	Trigger_CT[0]->SetThreshold(fpixelthres);
2207
2208	Trigger_CT[0]->Diskriminate();
2209
2210	//
2211	// look if in all the signals in the trigger signal branch
2212	// is a possible Trigger. Therefore we habe to diskriminate all
2213	// the simulated analog signals (Method Diskriminate in class
2214	// MTrigger). We look simultanously for the moments at which
2215	// there are more than TRIGGER_MULTI pixels above the
2216	// CHANNEL_THRESHOLD.
2217	//
2218
2219	// Set trigger flags to zero
2220	Lev0=0;
2221	Lev1=0;
2222
2223	// loop over multiplicity of trigger configuration
2224	for (imulticount = Trigger_loop_lmult;
2225	imulticount <= Trigger_loop_umult;
2226	imulticount++)
2227	{
2228	Trigger_CT[0]->SetMultiplicity(imulticount);
2229	Trigger_CT[0]->ClearZero();
2230
2231	Lev0=(Short_t) Trigger_CT[0]->ZeroLevel();
2232	if (Lev0>0 \|\| Write_All_Event_Headers \|\| btrigger)
2233	{
2234
2235	// loop over topologies
2236	for(itopocount=Trigger_loop_ltop;
2237	itopocount<=Trigger_loop_utop;
2238	itopocount++)
2239	{
2240	Lev1=0;
2241
2242	if(itopocount==0 && imulticount>7)
2243	continue;
2244
2245	//COBB if(itopocount==2 && imulticount<3) continue;
2246	// It only makes to look for a different topology
2247	// if there are 3 or more N pixels.
2248	if(imulticount<3)
2249	Trigger_CT[0]->SetTopology(1);
2250	else
2251	{
2252	// We should be careful that topologies are sort from
2253	// the less to the more restrictive one.
2254	Trigger_CT[0]->SetTopology(isorttopo[itopocount]);
2255	}
2256	Trigger_CT[0]->ClearFirst();
2257
2258	//
2259	// Start the First Level Trigger simulation
2260	//
2261	if(Lev0!=0)
2262	Lev1=Trigger_CT[0]->FirstLevel();
2263	if(Lev1>0) {
2264	btrigger= 1;
2265	ntriggerloop[ithrescount]
2266	[imulticount-Trigger_loop_lmult]
2267	[itopocount-Trigger_loop_ltop]++;
2268	}
2269
2270	Lev0=1;
2271	Int_t NumImages = Lev1;
2272	if(Lev1==0 && (Write_All_Event_Headers \|\| btrigger))
2273	{
2274	btrigger= 1;
2275	NumImages=1;
2276	Lev0=0;
2277	}
2278
2279	for (Int_t ii=0;ii<NumImages;ii++)
2280	{
2281	if (Write_McTrig)
2282	{
2283	McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
2284	McTrig[iconcount]->
2285	SetTime(Trigger_CT[0]->GetFirstLevelTime(ii),ii+1);
2286	Trigger_CT[0]->GetMapDiskriminator(trigger_map);
2287	McTrig[iconcount]->SetMapPixels(trigger_map,ii);
2288	}
2289	//
2290	// fill inside the class fadc the member output
2291	//
2292
2293	Fadc_CT[0]->TriggeredFadc(Trigger_CT[0]->
2294	GetFirstLevelTime(ii));
2295
2296	if( Write_RawEvt )
2297	{
2298	//
2299	// Fill the header of this event
2300	//
2301
2302	EvtHeader[iconcount]->
2303	FillHeader( (UInt_t) (ntshow + nshow),0);
2304
2305	// Set trigger pattern. WARNING: the trigger pattern of real
2306	// data is inverted for the MRawRunHeader.fFormatVersion=5
2307	// which is the one we now set in the camera outputs. Hence
2308	/// we use the binary inversion operator "~" below:
2309
2310	EvtHeader[iconcount]->
2311	SetTriggerPattern((UInt_t)~MTriggerPattern::kTriggerLvl1);
2312
2313	// fill pixel information
2314	if (Lev1){
2315	if (addElecNoise) Fadc_CT[0]->DigitalNoise();
2316	for(UInt_t i=0;
2317	i<((MGeomCam*)(camgeom.UncheckedAt(0)))
2318	->GetNumPixels();i++){
2319	//
2320	// AM 15 01 2004: commented out "continue"
2321	// statement, so that also pixels with no
2322	// C-photons will be written to the output
2323	// in case the camera is run with no noise.
2324	// if(!Fadc_CT[0]->IsPixelUsed(i)) continue;
2325
2326	for (j=0;j<FADC_slices_written;j++){
2327	fadcValues->AddAt(Fadc_CT[0]->
2328	GetFadcSignal(i,j),j);
2329	fadcValuesLow->AddAt(Fadc_CT[0]->
2330	GetFadcLowGainSignal(i,j),j);
2331	}
2332	EvtData[iconcount]->AddPixel(i,fadcValues,0);
2333	EvtData[iconcount]->AddPixel(i,fadcValuesLow,kTRUE);
2334	}
2335	}
2336	}
2337	}
2338	//
2339	// Increase counter of analised trigger conditions
2340	//
2341	iconcount++;
2342	}
2343	}
2344	else{
2345	break;
2346	}
2347	}
2348	if (!btrigger) break;
2349	}
2350	if (btrigger){
2351
2352	//
2353	// fill the MMcEvt with all information
2354	//
2355
2356	if(!flagstoring)
2357	nstoredevents++;
2358	flagstoring = 1;
2359
2360	if (Write_McEvt) {
2361	Float_t ftime, ltime;
2362	if (reflector_file_version<6){
2363	mcevth[0].get_times(&ftime, &ltime);
2364
2365	McEvtBasic[0]->Fill((MMcEvt::ParticleId_t) mcevth[0].get_primary(),
2366	mcevth[0].get_energy(), impactD[0],
2367	phiCT[0], thetaCT[0]);
2368
2369	McEvt[0]->Fill( 0,
2370	(UShort_t) mcevth[0].get_primary() ,
2371	mcevth[0].get_energy(),
2372	-1.0,
2373	-1.0,
2374	-1.0,
2375	mcevth[0].get_theta(),
2376	mcevth[0].get_phi(),
2377	mcevth[0].get_core(),
2378	coreX,
2379	coreY,
2380	impactD[0],
2381	phiCT[0],
2382	thetaCT[0],
2383	ftime,
2384	ltime,
2385	0,
2386	0,
2387	0,
2388	0,
2389	0,
2390	0,
2391	0,
2392	(UInt_t)mcevth[0].get_CORSIKA(),
2393	(UInt_t)mcevth[0].get_AtmAbs(),
2394	(UInt_t)(mcevth[0].get_MirrAbs()+
2395	mcevth[0].get_OutOfMirr()+
2396	mcevth[0].get_BlackSpot()),
2397	(UInt_t) ncph[0],
2398	(UInt_t) inumphe_CT[0],
2399	(UInt_t) inumphensb[0]+inumphe_CT[0],
2400	-1.0,
2401	-1.0,
2402	-1.0,
2403	fadc_jitter[0]);
2404	}
2405	else{
2406	Float_t Nmax, t0, tmax, a, b, c, chi2;
2407	mcevth_2[0].get_times(&ftime, &ltime);
2408	chi2=mcevth_2[0].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2409
2410	McEvtBasic[0]->Fill((MMcEvt::ParticleId_t) mcevth_2[0].get_primary(),
2411	mcevth_2[0].get_energy(), impactD[0],
2412	mcevth_2[0].get_phi_CT(),
2413	mcevth_2[0].get_theta_CT());
2414
2415	McEvt[0]->Fill((UInt_t) mcevth_2[0].get_evt_number(),
2416	(UShort_t) mcevth_2[0].get_primary() ,
2417	mcevth_2[0].get_energy(),
2418	mcevth_2[0].get_thick0(),
2419	mcevth_2[0].get_first_target(),
2420	mcevth_2[0].get_z_first_int(),
2421	mcevth_2[0].get_theta(),
2422	mcevth_2[0].get_phi(),
2423	mcevth_2[0].get_core(),
2424	coreX,
2425	coreY,
2426	impactD[0],
2427	mcevth_2[0].get_phi_CT(),
2428	mcevth_2[0].get_theta_CT(),
2429	ftime,
2430	ltime,
2431	Nmax,
2432	t0,
2433	tmax,
2434	a,
2435	b,
2436	c,
2437	chi2,
2438	(UInt_t)mcevth_2[0].get_CORSIKA(),
2439	(UInt_t)mcevth_2[0].get_AtmAbs(),
2440	(UInt_t)(mcevth_2[0].get_MirrAbs()+
2441	mcevth_2[0].get_OutOfMirr()+
2442	mcevth_2[0].get_BlackSpot()),
2443	(UInt_t) ncph[0],
2444	(UInt_t) inumphe_CT[0],
2445	(UInt_t) inumphensb[0]+inumphe_CT[0],
2446	mcevth_2[0].get_ElecFraction(),
2447	mcevth_2[0].get_MuonFraction(),
2448	mcevth_2[0].get_OtherFraction(),
2449	fadc_jitter[0]);
2450	}
2451	}
2452	// Fill the Tree with the current leaves of each branch
2453	i=EvtTree.Fill() ;
2454
2455	// Clear the branches
2456	if(Write_McTrig){
2457	for(int i=0;i<numberBranches;i++){
2458	McTrig[i]->Clear() ;
2459	}
2460	}
2461	if( Write_RawEvt ){
2462	for(int i=0;i<numberBranches;i++){
2463	EvtHeader[i]->Clear() ;
2464	EvtData[i]->ResetPixels (0, 0);
2465	}
2466	}
2467	if (Write_McEvt)
2468	{
2469	McEvt[0]->Clear();
2470	McEvtBasic[0]->Clear();
2471	}
2472	}
2473	}
2474
2475	// We study a single trigger condition
2476	else {
2477
2478	// Set to zero the flag to know if some conditon has triggered
2479	btrigger=0;
2480	flagstoring = 0;
2481
2482	for(int ict = 0; ict < ct_Number; ict++)
2483	{
2484
2485	// Setting trigger conditions
2486	Trigger_CT[ict]->SetMultiplicity(Trigger_multiplicity[ict]);
2487	Trigger_CT[ict]->SetTopology(Trigger_topology[ict]);
2488	for (int i=0;i<ct_NPixels;i++)
2489	fpixelthres[i]=qThreshold[ict][i];
2490
2491	// Rise the discrimnator threshold to avoid huge rates
2492	if(riseDiskThres>0.0 && simulateNSB && inumphe >= nphe2NSB)
2493	for(int ii=0;ii<ct_NPixels;ii++)
2494	{
2495	if( nsb_phepns_rotated[ict][ii]>riseDiskThres)
2496	fpixelthres[ii]=secureDiskThres;
2497	}
2498
2499	Trigger_CT[ict]->SetThreshold(fpixelthres);
2500
2501	Trigger_CT[ict]->Diskriminate() ;
2502
2503	//
2504	// Look if in all the signals in the trigger signal branch
2505	// is a possible Trigger. Therefore we have to discriminate all
2506	// the simulated analog signals (Method Diskriminate in class
2507	// MTrigger). We look simultaneously for the moments at which
2508	// there are more than TRIGGER_MULTI pixels above the
2509	// CHANNEL_THRESHOLD.
2510	//
2511
2512	Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ;
2513
2514	Lev1MT[ict] = 0 ;
2515
2516	//
2517	// Start the First Level Trigger simulation
2518	//
2519
2520	if ( Lev0MT[ict] > 0 /\|\| Write_All_Event_Headers/)
2521	Lev1MT[ict]= Trigger_CT[ict]->FirstLevel();
2522
2523	if (Lev1MT[ict]>0)
2524	++ntrigger[ict];
2525
2526	}
2527
2528	Int_t NumImages = 0;
2529	Int_t CT_triggered=0;
2530	for(int ict=0;ict<ct_Number;ict++)
2531	{
2532	if(NumImages==0 && Lev1MT[ict]>0)
2533	CT_triggered=ict;
2534	NumImages = (NumImages>=Lev1MT[ict]) ? NumImages : 1;
2535	Lev0MT[ict]=1;
2536
2537	if (Lev1MT[ict]==0 && Write_All_Event_Headers)
2538	{
2539	NumImages=1;
2540	Lev0MT[ict]=0;
2541	}
2542	}
2543
2544	for(int ict=0;ict<ct_Number;ict++){
2545	for(Int_t ii=0;ii<NumImages;ii++){
2546
2547	btrigger=1;
2548
2549	// Loop over different level one triggers
2550
2551	//
2552	// fill inside class fadc the member output
2553	//
2554	if(Lev1MT[ict]>0)
2555	Fadc_CT[ict]->
2556	TriggeredFadc(Trigger_CT[ict]->GetFirstLevelTime(ii));
2557	else
2558	Fadc_CT[ict]->
2559	TriggeredFadc(Trigger_CT[CT_triggered]->GetFirstLevelTime(ii));
2560
2561	if(!flagstoring)
2562	nstoredevents++;
2563	flagstoring = 1;
2564
2565	if (Write_McTrig){
2566	McTrig[ict]->SetFirstLevel (Lev1MT[ict]);
2567	McTrig[ict]
2568	->SetTime(Trigger_CT[ict]->GetFirstLevelTime(ii),ii+1);
2569	Trigger_CT[ict]->GetMapDiskriminator(trigger_map);
2570	McTrig[ict]->SetMapPixels(trigger_map,ii);
2571	}
2572
2573	// Fill Evt information
2574
2575	if (Write_RawEvt){
2576
2577	//
2578	// Fill the header of this event
2579	//
2580
2581	EvtHeader[ict]
2582	->FillHeader ( (UInt_t) (ntshow + nshow) , 0 ) ;
2583
2584
2585	// Set trigger pattern. WARNING: the trigger pattern of real
2586	// data is inverted for the MRawRunHeader.fFormatVersion=5
2587	// which is the one we now set in the camera outputs. Hence
2588	/// we use the binary inversion operator "~" below:
2589	EvtHeader[ict]->SetTriggerPattern((UInt_t)~MTriggerPattern::kTriggerLvl1);
2590
2591	// Fill pixel information
2592	// AM 17/2/2005: added condition on inumphe. Noise is not generated for
2593	// events with less than phe2NSB photoelectrons, and then it is better not
2594	// to write the images of those events to the output, so that we avoid them
2595	// to be processed (with no noise) later in the chain.
2596
2597	if (Lev1MT[ict] && inumphe >= nphe2NSB)
2598	{
2599	if (addElecNoise)
2600	Fadc_CT[ict]->DigitalNoise();
2601	for(UInt_t i=0;
2602	i<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
2603	i++){
2604	//
2605	// AM 15 01 2004 Commented out "continue" statement, so that also pixels
2606	// containing no C-photons will be written to the output in case of running
2607	// camera with no noise added to the signal.
2608	// if(!Fadc_CT[ict]->IsPixelUsed(i)) continue;
2609	//
2610	for (j = 0; j < FADC_slices_written; j++)
2611	{
2612	fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(i,j),j);
2613	fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(i,j),j);
2614	}
2615	EvtData[ict]->AddPixel(i,fadcValues,0);
2616	EvtData[ict]->AddPixel(i,fadcValuesLow,kTRUE);
2617	}
2618	}
2619	}
2620	}
2621	//
2622	// if a first level trigger occurred, then
2623	// 1. do some other stuff (not implemented)
2624	// 2. start the gui tool
2625
2626	if(FADC_Scan){
2627	if ( Lev0MT[ict] > 0 ) {
2628	Fadc_CT[ict]->ShowSignal( McEvt[ict], (Float_t) 60. ) ;
2629	}
2630	}
2631
2632	if(Trigger_Scan){
2633	if ( Lev0MT[ict] > 0 ) {
2634	Trigger_CT[ict]->ShowSignal(McEvt[ict]) ;
2635	}
2636	}
2637
2638	} // end CT loop
2639
2640	// If there is trigger in some telescope or we store all showers
2641	if(btrigger){
2642	if (Write_McEvt){
2643	//
2644	// fill the MMcEvt with all information
2645	//
2646
2647	for (int ict=0;ict<ct_Number;ict++){
2648	Float_t ftime, ltime;
2649	if (reflector_file_version<6){
2650	mcevth[ict].get_times(&ftime, &ltime);
2651
2652	McEvtBasic[ict]->Fill((MMcEvt::ParticleId_t) mcevth[0].get_primary(),
2653	mcevth[ict].get_energy(), impactD[ict],
2654	phiCT[ict], thetaCT[ict]);
2655
2656	McEvt[ict]->Fill( 0,
2657	(UShort_t) mcevth[0].get_primary() ,
2658	mcevth[ict].get_energy(),
2659	-1.0,
2660	-1.0,
2661	-1.0,
2662	mcevth[ict].get_theta(),
2663	mcevth[ict].get_phi(),
2664	mcevth[ict].get_core(),
2665	coreX,
2666	coreY,
2667	impactD[ict],
2668	phiCT[ict],
2669	thetaCT[ict],
2670	ftime,
2671	ltime,
2672	0,
2673	0,
2674	0,
2675	0,
2676	0,
2677	0,
2678	0,
2679	(UInt_t)mcevth[ict].get_CORSIKA(),
2680	(UInt_t)mcevth[ict].get_AtmAbs(),
2681	(UInt_t)(mcevth[ict].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()),
2682	(UInt_t) ncph[ict],
2683	(UInt_t) inumphe_CT[ict],
2684	(UInt_t) inumphensb[ict]+inumphe_CT[ict],
2685	-1.0,
2686	-1.0,
2687	-1.0,
2688	fadc_jitter[ict]);
2689	}
2690	else{
2691	Float_t Nmax, t0, tmax, a, b, c, chi2;
2692	mcevth_2[ict].get_times(&ftime, &ltime);
2693	chi2=mcevth_2[ict].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2694
2695	McEvtBasic[ict]->Fill((MMcEvt::ParticleId_t)mcevth_2[ict].get_primary(),
2696	mcevth_2[ict].get_energy(), impactD[ict],
2697	mcevth_2[ict].get_phi_CT(),
2698	mcevth_2[ict].get_theta_CT());
2699
2700
2701	McEvt[ict]->Fill( (UInt_t) mcevth_2[ict].get_evt_number(),
2702	(UShort_t) mcevth_2[ict].get_primary() ,
2703	mcevth_2[ict].get_energy(),
2704	mcevth_2[ict].get_thick0(),
2705	mcevth_2[ict].get_first_target(),
2706	mcevth_2[ict].get_z_first_int(),
2707	mcevth_2[ict].get_theta(),
2708	mcevth_2[ict].get_phi(),
2709	mcevth_2[ict].get_core(),
2710	coreX,
2711	coreY,
2712	impactD[ict],
2713	mcevth_2[ict].get_phi_CT(),
2714	mcevth_2[ict].get_theta_CT(),
2715	ftime,
2716	ltime,
2717	Nmax,
2718	t0,
2719	tmax,
2720	a,
2721	b,
2722	c,
2723	chi2,
2724	(UInt_t)mcevth_2[ict].get_CORSIKA(),
2725	(UInt_t)mcevth_2[ict].get_AtmAbs(),
2726	(UInt_t) (mcevth_2[ict].get_MirrAbs()+mcevth_2[ict].get_OutOfMirr()+mcevth_2[ict].get_BlackSpot()),
2727	(UInt_t) ncph[ict],
2728	(UInt_t) inumphe_CT[ict],
2729	(UInt_t) inumphensb[ict]+inumphe_CT[ict],
2730	mcevth_2[ict].get_ElecFraction(),
2731	mcevth_2[ict].get_MuonFraction(),
2732	mcevth_2[ict].get_OtherFraction(),
2733	fadc_jitter[ict]);
2734	}
2735	}
2736	}
2737	// We do not count photons out of the camera.
2738
2739
2740	//
2741	// write it out to the file outfile
2742	//
2743
2744	EvtTree.Fill() ;
2745	}
2746
2747	// clear all
2748	for(int ict=0;ict<ct_Number;ict++){
2749	if (Write_RawEvt) EvtHeader[ict]->Clear() ;
2750	if (Write_RawEvt) EvtData[ict]->ResetPixels(0,0);
2751	if (Write_McTrig) McTrig[ict]->Clear() ;
2752	if (Write_McEvt)
2753	{
2754	McEvt[ict]->Clear() ;
2755	McEvtBasic[ict]->Clear();
2756	}
2757	}
2758	}
2759
2760	#ifdef __DEBUG__
2761	printf("\n");
2762
2763	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
2764
2765	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
2766
2767	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
2768
2769	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
2770
2771	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
2772	(int)pixels[ici][icj][PIXNUM],
2773	pixels[ici][icj][PIXX],
2774	pixels[ici][icj][PIXY],
2775	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
2776
2777	}
2778	}
2779	}
2780	}
2781
2782	for (int i=0;
2783	i<((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels(); ++i) {
2784	printf("%d (%d): ", i, npixneig[i]);
2785	for (j=0; j<npixneig[i]; ++i)
2786	printf(" %d", pixneig[i][j]);
2787	printf("\n");
2788	}
2789
2790	#endif // __DEBUG__
2791
2792
2793	// We search the maximum impact parameter fo the simualted showers
2794	maxpimpact=maxpimpact<impactD[0]?impactD[0]:maxpimpact;
2795
2796	// look for the next event
2797
2798	for(int ict=0;ict<ct_Number;ict++)
2799	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2800
2801	} // end while there is a next event
2802
2803	if( !isA( flag, FLAG_END_OF_RUN )){
2804	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
2805	break;
2806	}
2807	else { // found end of run
2808	ntshow += nshow;
2809	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
2810
2811	//fread( flag, SIZE_OF_FLAGS, 1, inputfile );
2812	for(int ict=0;ict<ct_Number;ict++)
2813	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2814
2815	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
2816	log(SIGNATURE, "End of file . . .\n");
2817	still_in_loop = FALSE;
2818	log(SIGNATURE, "Reading ASCII files at the end of the reflector file. . .\n");
2819	for(int ict=0;ict<ct_Number;ict++){
2820	read_ascii(inputfile[ict], McConfigRunHeader[ict]);
2821	McConfigRunHeader[ict]->SetMissPointingX(missP_x);
2822	McConfigRunHeader[ict]->SetMissPointingY(missP_y);
2823
2824	McConfigRunHeader[ict]->SetMirrorFraction(mirror_frac[ict]);
2825
2826	if ( Spotsigma > 0.)
2827	{
2828	Float_t ref_spotsigma = McConfigRunHeader[ict]->GetPointSpread();
2829	Float_t newsigma = sqrt(ref_spotsigma * ref_spotsigma +
2830	Spot_x* Spot_x);
2831	McConfigRunHeader[ict]->SetPointSpreadX(newsigma);
2832	newsigma = sqrt(ref_spotsigma * ref_spotsigma + Spot_y* Spot_y);
2833	McConfigRunHeader[ict]->SetPointSpreadY(newsigma);
2834	}
2835	}
2836	if ((! Data_From_STDIN) && ( !feof(inputfile[0]) )){
2837
2838	// we have concatenated input files.
2839	// get signature of the next part and check it.
2840
2841	if((reflector_file_version=check_reflector_file( inputfile[0] ))==FALSE){
2842	log(SIGNATURE, "Next file is not recognised as a reflector file.\n");
2843	log(SIGNATURE, "Stopping ...\n");
2844	break;
2845	}
2846
2847	}
2848
2849	for(int ict=0;ict<ct_Number;ict++)
2850	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2851
2852	} // end if found end of file
2853	} // end if found end of run
2854
2855	} // end if else found start of run
2856	} // end big while loop
2857
2858	//<@ Finally we should fill the McRunHeader
2859
2860	Float_t heights[10];
2861	time_t ltime;
2862	Float_t ftime;
2863	Float_t rnum;
2864	Float_t viewcone[2]={0,0};
2865
2866	get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
2867	if (reflector_file_version<6){
2868	mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2869	mcevth[0].get_phi_range(&shphimin,&shphimax);
2870	mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2871	mcevth[0].get_phi_range(&shphimin,&shphimax);
2872	corsika=UInt_t(mcevth[0].get_VersionPGM()*1000);
2873	for (int i=0; i< 10;i++)
2874	heights[i]=mcevth[0].get_HeightLev (i);
2875	rnum=mcevth[0].get_RunNumber();
2876	}
2877	else{
2878	mcevth_2[0].get_theta_range(&shthetamin, &shthetamax);
2879	mcevth_2[0].get_phi_range(&shphimin,&shphimax);
2880	corsika=UInt_t(mcevth_2[0].get_VersionPGM()*1000);
2881	for (int i=0; i< 10;i++)
2882	heights[i]=mcevth_2[0].get_HeightLev (i);
2883	rnum=mcevth_2[0].get_RunNumber();
2884	mcevth_2[0].get_viewcone(&viewcone[0],&viewcone[1]);
2885	}
2886
2887	if(!Trigger_Loop) icontrigger=0;
2888	time (&ltime);
2889	ftime = ((Float_t)ltime)/1000;
2890
2891	if (reflector_file_version<6)
2892	McRunHeader->Fill(rnum,
2893	(UInt_t) 0,
2894	mcevth[0].get_DateRun(),
2895	ftime,
2896	icontrigger,
2897	!Write_All_Event_Headers,
2898	Write_McEvt,
2899	Write_McTrig,
2900	Write_McFADC,
2901	Write_RawEvt,
2902	addElecNoise,
2903	ct_NPixels,
2904	(UInt_t)ntshow,
2905	(UInt_t)nstoredevents,
2906	0,
2907	sfRaH,
2908	sfRaM,
2909	sfRaS,
2910	sfDeD,
2911	sfDeM,
2912	sfDeS,
2913	meanNSB,
2914	shthetamax,
2915	shthetamin,
2916	shphimax,
2917	shphimin,
2918	maxpimpact,
2919	mcevth[0].get_CWaveLower(),
2920	mcevth[0].get_CWaveUpper(),
2921	mcevth[0].get_slope(),
2922	1,
2923	heights,
2924	corsika,
2925	(UInt_t)(reflector_file_version*100),
2926	(UInt_t)(VERSION*100),
2927	0);
2928	else
2929	McRunHeader->Fill(rnum,
2930	(UInt_t) 0,
2931	mcevth_2[0].get_DateRun(),
2932	ftime,
2933	icontrigger,
2934	!Write_All_Event_Headers,
2935	Write_McEvt,
2936	Write_McTrig,
2937	Write_McFADC,
2938	Write_RawEvt,
2939	addElecNoise,
2940	ct_NPixels,
2941	(UInt_t)ntshow,
2942	(UInt_t)nstoredevents,
2943	0,
2944	sfRaH,
2945	sfRaM,
2946	sfRaS,
2947	sfDeD,
2948	sfDeM,
2949	sfDeS,
2950	meanNSB,
2951	shthetamax,
2952	shthetamin,
2953	shphimax,
2954	shphimin,
2955	maxpimpact,
2956	mcevth_2[0].get_CWaveLower(),
2957	mcevth_2[0].get_CWaveUpper(),
2958	mcevth_2[0].get_slope(),
2959	1,
2960	heights,
2961	corsika,
2962	(UInt_t)(reflector_file_version*100),
2963	(UInt_t)(VERSION*100),
2964	0);
2965	// Fill some missing values for MRawRunHeader
2966
2967	RunHeader->SetRunNumber((UInt_t)rnum);
2968	RunHeader->SetNumEvents(nstoredevents);
2969
2970	// Fill MMcCorsikaRunHeader
2971
2972	Float_t constantC[50];
2973	mcrunh.get_constantC(&constantC[0]);
2974	Float_t constantCKA[40];
2975	mcrunh.get_constantCKA(&constantCKA[0]);
2976	Float_t constantCETA[5];
2977	mcrunh.get_constantCETA(&constantCETA[0]);
2978	Float_t constantCSTRBA[11];
2979	mcrunh.get_constantCSTRBA(&constantCSTRBA[0]);
2980	Float_t constantAATM[5];
2981	mcrunh.get_constantAATM(&constantAATM[0]);
2982	Float_t constantBATM[5];
2983	mcrunh.get_constantBATM(&constantBATM[0]);
2984	Float_t constantCATM[5];
2985	mcrunh.get_constantCATM(&constantCATM[0]);
2986	Float_t constantNFL[4];
2987	mcrunh.get_constantNFL(&constantNFL[0]);
2988
2989	if(reflector_file_version>5)
2990	McCorsikaRunHeader->Fill(rnum,
2991	mcrunh.get_date(),
2992	corsika,
2993	1,
2994	heights,
2995	mcevth_2[0].get_slope(),
2996	mcrunh.get_ELow(),
2997	mcrunh.get_EUpp(),
2998	mcrunh.get_EGS4(),
2999	mcrunh.get_NKG(),
3000	mcrunh.get_Ecutoffh(),
3001	mcrunh.get_Ecutoffm(),
3002	mcrunh.get_Ecutoffe(),
3003	mcrunh.get_Ecutoffg(),
3004	constantC,
3005	constantCKA,
3006	constantCETA,
3007	constantCSTRBA,
3008	constantAATM,
3009	constantBATM,
3010	constantCATM,
3011	constantNFL,
3012	viewcone,
3013	mcrunh.get_wobble(),
3014	mcrunh.get_atmophere()
3015	);
3016
3017	// Store qe for each PMT in output file
3018	TArrayF qe_pmt;
3019	TArrayF wav_pmt;
3020
3021	for(int ict=0;ict<ct_Number;ict++){
3022	McConfigRunHeader[ict]->InitSizePMTs(ct_NPixels);
3023	for(int i=0; i<(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();i++){
3024	McConfigRunHeader[ict]->SetPmtTimeJitter(pmt_jitter);
3025	McConfigRunHeader[ict]->AddPMT(i);
3026	MGeomPMT &pmt = McConfigRunHeader[ict]->GetPMT(i);
3027	qe_pmt.Set(pointsQE[ict],QE[ict][i][1]);
3028	wav_pmt.Set(pointsQE[ict],QElambda);
3029	pmt.SetArraySize(pointsQE[ict]);
3030	pmt.SetPMTContent(i,wav_pmt,qe_pmt);
3031	}
3032
3033	// Store Light Collection factors in the output file
3034	TArrayF theta_lc;
3035	TArrayF factor_lc;
3036	TArrayF factor_lc_outer;
3037
3038	theta_lc.Set(pointsWC,WC[0]);
3039	factor_lc.Set(pointsWC,WC[1]);
3040	factor_lc_outer.Set(pointsWC,WC_outer[1]);
3041
3042	McConfigRunHeader[ict]->SetLightCollection(theta_lc, factor_lc,
3043	factor_lc_outer);
3044
3045	}
3046
3047	// Fill the Header Tree with the current leaves of each branch
3048	HeaderTree.Fill() ;
3049
3050	//++
3051	// put the Event to the root file
3052	//--
3053
3054	outfile_temp.Write() ;
3055	outfile_temp.Close() ;
3056
3057	// close input file
3058
3059	if (Trigger_Loop){
3060	log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
3061	ntshow, ntcph[0] );
3062	datafile<<ntshow<<" event(s), with a total of "<<ntcph[0]<<" C.photons"<<endl;
3063	log( SIGNATURE, "Trigger Mode. \n");
3064	log( SIGNATURE, "Fraction of triggers: \n");
3065	datafile<<"Fraction of triggers: "<<endl;
3066	for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
3067	for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
3068	for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
3069	log( SIGNATURE, "Thres %5.1f, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
3070	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);
3071	datafile<<"Thres "<<fthrescount<<", Multi "<<imulticount<<", Topo"<<isorttopo[itopocount]<<": ";
3072	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;
3073	}
3074	}
3075	}
3076	}
3077	else{
3078	for(int ict=0;ict<ct_Number;ict++){
3079	log( SIGNATURE,
3080	"%d event(s), with a total of %d C.photons in CT %i (%s)\n",
3081	ntshow, ntcph[ict],ict,GeometryName[ict] );
3082	datafile<<ntshow<<" event(s), with a total of "<<ntcph[ict]
3083	<<" C.photons in CT "<<ict<<" ("<<GeometryName[ict]<<")"<<endl;
3084	log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
3085	((float)ntrigger[ict]) / ((float)ntshow) * 100.0, ntrigger[ict], ntshow);
3086	datafile<<"Fraction of triggers: "<<((float)ntrigger[ict]) / ((float)ntshow) * 100.0<<" ("<<ntrigger[ict]<<" out of "<<ntshow<<" )"<<endl;
3087	}
3088	}
3089
3090	// close files
3091
3092	log( SIGNATURE, "Closing files\n" );
3093
3094	if( ! Data_From_STDIN ){
3095	for(int ict=0;ict<ct_Number;ict++)
3096	fclose( inputfile[ict] );
3097	}
3098	datafile.close();
3099
3100	// program finished
3101
3102	log( SIGNATURE, "Done.\n");
3103
3104	return( 0 );
3105	}
3106	//!@}
3107
3108	// @T \newpage
3109
3110	//!@subsection Functions definition.
3111
3112	//!-----------------------------------------------------------
3113	// @name present
3114	//
3115	// @desc Make some presentation
3116	//
3117	// @date Sat Jun 27 05:58:56 MET DST 1998
3118	//------------------------------------------------------------
3119	// @function
3120
3121	//!@{
3122	void
3123	present(void)
3124	{
3125	cout << "##################################################\n"
3126	<< SIGNATURE << '\n' << '\n'
3127	<< "Processor of the reflector output\n"
3128	<< "J. C. Gonzalez, Jun 1998\n"
3129	<< "O. Blanch, A. Moralejo, 2004\n"
3130	<< "##################################################\n\n"
3131	<< flush ;
3132	}
3133	//!@}
3134
3135
3136	//!-----------------------------------------------------------
3137	// @name usage
3138	//
3139	// @desc show help
3140	//
3141	// @date Tue Dec 15 16:23:30 MET 1998
3142	//------------------------------------------------------------
3143	// @function
3144
3145	//!@{
3146	void
3147	usage(void)
3148	{
3149	present();
3150	cout << "\nusage ::\n\n"
3151	<< "\t camera "
3152	<< " [ -@ paramfile ] "
3153	<< " [ -h ] "
3154	<< "\n\n or \n\n"
3155	<< "\t camera < paramfile"
3156	<< "\n\n";
3157	exit(0);
3158	}
3159	//!@}
3160
3161
3162	//!-----------------------------------------------------------
3163	// @name log
3164	//
3165	// @desc function to send log information
3166	//
3167	// @var funct Name of the caller function
3168	// @var fmt Format to be used (message)
3169	// @var ... Other information to be shown
3170	//
3171	// @date Sat Jun 27 05:58:56 MET DST 1998
3172	//------------------------------------------------------------
3173	// @function
3174
3175	//!@{
3176	void
3177	log(const char funct, char fmt, ...)
3178	{
3179	va_list args;
3180
3181	// Display the name of the function that called error
3182	printf("[%s]: ", funct);
3183
3184	// Display the remainder of the message
3185	va_start(args, fmt);
3186	vprintf(fmt, args);
3187	va_end(args);
3188	}
3189	//!@}
3190
3191
3192	//!-----------------------------------------------------------
3193	// @name error
3194	//
3195	// @desc function to send an error message, and abort the program
3196	//
3197	// @var funct Name of the caller function
3198	// @var fmt Format to be used (message)
3199	// @var ... Other information to be shown
3200	//
3201	// @date Sat Jun 27 05:58:56 MET DST 1998
3202	//------------------------------------------------------------
3203	// @function
3204
3205	//!@{
3206	void
3207	error(const char funct, char fmt, ...)
3208	{
3209	va_list args;
3210
3211	// Display the name of the function that called error
3212	fprintf(stdout, "ERROR in %s: ", funct);
3213
3214	// Display the remainder of the message
3215	va_start(args, fmt);
3216	vfprintf(stdout, fmt, args);
3217	va_end(args);
3218
3219	perror(funct);
3220
3221	exit(1);
3222	}
3223	//!@}
3224
3225
3226	//!-----------------------------------------------------------
3227	// @name isA
3228	//
3229	// @desc returns TRUE(FALSE), if the flag is(is not) the given
3230	//
3231	// @var s1 String to be searched
3232	// @var flag Flag to compare with string s1
3233	// @return TRUE: both strings match; FALSE: oth.
3234	//
3235	// @date Wed Jul 8 15:25:39 MET DST 1998
3236	//------------------------------------------------------------
3237	// @function
3238
3239	//!@{
3240	int
3241	isA( char * s1, const char * flag ) {
3242	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
3243	}
3244	//!@}
3245
3246
3247	//!-----------------------------------------------------------
3248	// @name read_QE
3249	//
3250	// @desc read QE data
3251	//
3252	// @date thu 5 17:59:57 CEST 2002
3253	//------------------------------------------------------------
3254	// @function
3255
3256	//!@{
3257	void
3258	read_QE(char fname[256], int ict){
3259	ifstream qefile;
3260	char line[LINE_MAX_LENGTH];
3261	int i, j, icount;
3262	float qe;
3263
3264	//------------------------------------------------------------
3265	// second, pixels' QE
3266
3267	// try to open the file
3268
3269	log("read_QE", "Opening the file \"%s\" . . .\n", fname);
3270
3271	qefile.open( fname );
3272
3273	// if it is wrong or does not exist, exit
3274
3275	if ( qefile.bad() )
3276	error( "read_QE", "Cannot open \"%s\". Exiting.\n", fname );
3277
3278	// read file
3279
3280	log("read_QE", "Reading data . . .\n");
3281
3282	i=-1;
3283	icount = 0;
3284
3285	while ( ! qefile.eof() ) {
3286
3287	// get line from the file
3288
3289	qefile.getline(line, LINE_MAX_LENGTH);
3290
3291	// skip if comment
3292
3293	if ( *line == '#' )
3294	continue;
3295
3296	// if it is the first valid value, it is the number of QE data points
3297
3298	if ( i < 0 ) {
3299
3300	// get the number of datapoints
3301
3302	sscanf(line, "%d", &pointsQE[ict]);
3303
3304	// allocate memory for the table of QEs
3305
3306	QE[ict] = new float ** [ct_NPixels];
3307
3308	for ( i=0; i<ct_NPixels; ++i ) {
3309	QE[ict][i] = new float * [2];
3310	QE[ict][i][0] = new float[pointsQE[ict]];
3311	QE[ict][i][1] = new float[pointsQE[ict]];
3312	}
3313
3314	QElambda = new float [pointsQE[ict]];
3315
3316	for ( i=0; i<pointsQE[ict]; ++i ) {
3317	qefile.getline(line, LINE_MAX_LENGTH);
3318	sscanf(line, "%f", &QElambda[i]);
3319	}
3320
3321	i=0;
3322
3323	continue;
3324	}
3325
3326	// get the values (num-pixel, num-datapoint, QE-value)
3327
3328	if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
3329	break;
3330
3331	if ( (i < ct_NPixels) && (i > -1) &&
3332	((j-1) < pointsQE[ict]) && ((j-1) > -1) ) {
3333	QE[ict][i][0][j-1] = QElambda[j-1];
3334	QE[ict][i][1][j-1] = qe;
3335	}
3336
3337	if ( i > ct_NPixels) break;
3338
3339	icount++;
3340
3341	}
3342
3343	if(icount/pointsQE[ict] < ct_NPixels){
3344	error( "read_QE", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
3345	icount/pointsQE[ict], ct_NPixels );
3346	}
3347
3348	// close file
3349
3350	qefile.close();
3351
3352	// test QE
3353
3354	for(icount=0; icount< ct_NPixels; icount++){
3355	for(i=0; i<pointsQE[ict]; i++){
3356	if( QE[ict][icount][0][i] < 100. \|\| QE[ict][icount][0][i] > 1000. \|\|
3357	QE[ict][icount][1][i] < 0. \|\| QE[ict][icount][1][i] > 100.){
3358	error( "read_QE", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
3359	icount, i, QE[ict][icount][0][i], QE[ict][icount][1][i] );
3360	}
3361	}
3362	}
3363
3364	// end
3365
3366	log("read_QE", "Done.\n");
3367	}
3368	//!@}
3369
3370	//!-----------------------------------------------------------
3371	// @name read_WC
3372	//
3373	// @desc read WC data
3374	//
3375	// @date thu 5 17:59:57 CEST 2002
3376	//------------------------------------------------------------
3377	// @function
3378
3379	//!@{
3380	void
3381	read_WC(void){
3382	ifstream wcfile;
3383	char line[LINE_MAX_LENGTH];
3384	int i;
3385
3386	//------------------------------------------------------------
3387	// Read Light Collection data
3388
3389	// try to open the file
3390
3391	log("read_WC", "Opening the file \"%s\" . . .\n", WC_FILE);
3392
3393	wcfile.open( WC_FILE );
3394
3395	// if it is wrong or does not exist, exit
3396
3397	if ( wcfile.bad() )
3398	error( "read_WC", "Cannot open \"%s\". Exiting.\n", WC_FILE );
3399
3400	// read file
3401
3402	log("read_WC", "Reading data . . .\n");
3403
3404	// get line from the file
3405
3406	do
3407	wcfile.getline(line, LINE_MAX_LENGTH);
3408	while (line[0] == '#');
3409
3410	// get the number of datapoints
3411
3412	sscanf(line, "%d", &pointsWC);
3413
3414	// allocate memory for the table of QEs
3415
3416	WC = new float * [2];
3417	WC[0] = new float[pointsWC];
3418	WC[1] = new float[pointsWC];
3419
3420	for ( i=0; i<pointsWC; ++i ) {
3421	wcfile.getline(line, LINE_MAX_LENGTH);
3422	sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
3423	}
3424
3425	//
3426	// Now read info for outer pixels
3427	//
3428	WC_outer = new float * [2];
3429	WC_outer[0] = new float[pointsWC];
3430	WC_outer[1] = new float[pointsWC];
3431
3432	do
3433	wcfile.getline(line, LINE_MAX_LENGTH);
3434	while (line[0] == '#');
3435
3436	if (wcfile.eof())
3437	{
3438	log("read_WC", "ERROR. Missing data for outer pixels in file \"%s\"...\n",WC_FILE);
3439	log("read_WC", "EXITING camera\n");
3440	exit(-1);
3441	}
3442
3443	sscanf(line, "%f %f", &WC_outer[0][0], &WC_outer[1][0]);
3444	for ( i=1; i<pointsWC; ++i ) {
3445	wcfile.getline(line, LINE_MAX_LENGTH);
3446	sscanf(line, "%f %f", &WC_outer[0][i], &WC_outer[1][i]);
3447	}
3448
3449	// close file
3450
3451	wcfile.close();
3452
3453	// read
3454
3455	log("read_WC", "Done.\n");
3456	}
3457	//!@}
3458
3459
3460	//!-----------------------------------------------------------
3461	// @name read_ascii
3462	//
3463	// @desc read ascii configuration files used by the reflector
3464	//
3465	// @date tue dec 10 17:14:10 CET 2002
3466	//------------------------------------------------------------
3467	// @function
3468
3469	//!@{
3470	void
3471	read_ascii(FILE sp, MMcConfigRunHeader config)
3472	{
3473	Float_t radius = -1.0;
3474	Float_t focal = -1.0;
3475	Float_t point = -1.0;
3476	Float_t spot = -1.0;
3477	Float_t camwidth = -1.0;
3478
3479	Int_t imir;
3480	Float_t f,sx,sy,x,y,z,thetan,phin,xn,yn,zn;
3481	Float_t dx,dy;
3482	Int_t nummir, numref;
3483	Float_t wav,ref;
3484
3485	Char_t token[40];
3486	Char_t line[511];
3487	Char_t flag;
3488
3489	while(1){
3490	if((flag=fgetc(sp))==EOF)
3491	break;
3492	if (flag == '\n') // skip empty lines
3493	continue;
3494
3495	fgets(&line[1],500,sp);
3496	line[0]=flag;
3497
3498	if ( strstr(line, "# reflectivity file") == line ) {
3499	while (strstr(line, "# number of datapoints") != line)
3500	fgets(line,500,sp);
3501
3502	fgets(line,500,sp);
3503	sscanf(line,"%i",&numref);
3504
3505	TArrayF wavarray(numref);
3506	TArrayF refarray(numref);
3507
3508	while (strstr(line, "# datapoints") != line)
3509	fgets(line,500,sp);
3510
3511	for(int i=0; i<numref;i++){
3512	fgets(line,500,sp);
3513	if (line[0] == '#')
3514	{
3515	i--;
3516	continue;
3517	}
3518	sscanf(line,"%f %f",&wav,&ref);
3519	wavarray[i]=wav;
3520	refarray[i]=ref;
3521	}
3522	for (int j=0; j<nummir;j++){
3523
3524	MGeomMirror &mirror = config->GetMirror(j);
3525	mirror.SetArraySize(numref);
3526	mirror.SetReflectivity(wavarray, refarray);
3527	}
3528	continue;
3529	}
3530	if (line[0]== '#')
3531	continue;
3532	if (strstr(line, "type") == line)
3533	continue;
3534	if (strstr(line, "focal_distance") == line){
3535	sscanf(line,"%s %f",token,&focal);
3536	continue;
3537	}
3538	if (strstr(line, "point_spread") == line){
3539	sscanf(line,"%s %f",token,&point);
3540	continue;
3541	}
3542	if (strstr(line, "black_spot") == line){
3543	sscanf(line,"%s %f",token,&spot);
3544	continue;
3545	}
3546	if (strstr(line, "camera_width") == line){
3547	sscanf(line,"%s %f",token,&camwidth);
3548	continue;
3549	}
3550	//
3551	// Skip obsolete magic.def entries:
3552	//
3553	if (strstr(line, "n_pixels") == line)
3554	continue;
3555	if (strstr(line, "pixel_width") == line)
3556	continue;
3557	if (strstr(line, "n_centralpixels") == line)
3558	continue;
3559	if (strstr(line, "n_gappixels") == line)
3560	continue;
3561
3562	if (strstr(line, "n_mirrors") == line){
3563	sscanf(line,"%s %i",token,&nummir);
3564	config->InitSizeMirror(nummir);
3565	continue;
3566	}
3567	if (strstr(line, "r_mirror") == line){
3568	sscanf(line,"%s %f",token,&radius);
3569	continue;
3570	}
3571	if (strstr(line, "define_mirrors") == line){
3572	for(int i=0;i<nummir;i++){
3573	fgets(line,500,sp);
3574	sscanf(line,"%i %f %f %f %f %f %f %f %f %f %f %f",
3575	&imir,&f,&sx,&sy,&x,&y,&z,&thetan,&phin,&xn,&yn,&zn);
3576	config->AddMirror(i);
3577	MGeomMirror &mirror = config->GetMirror(i);
3578	mirror.SetMirrorContent(imir,f,sx,sy,x,y,z,thetan,phin,xn,yn,zn);
3579	}
3580	fgets(line,500,sp);
3581
3582	while ( ! strstr(line, "axis deviation"))
3583	fgets(line,500,sp);
3584
3585	for(int i=0;i<nummir;i++){
3586	fgets(line,500,sp);
3587	if (line[0] == '#')
3588	{
3589	i--;
3590	continue;
3591	}
3592	sscanf(line,"%f %f",&dx,&dy);
3593	MGeomMirror &mirror = config->GetMirror(i);
3594	mirror.SetMirrorDeviations(dx,dy);
3595	}
3596	continue;
3597	}
3598	}
3599	config->SetMagicDef(radius, focal, point, spot, camwidth);
3600	}
3601
3602
3603	//!-----------------------------------------------------------
3604	// @name igen_pixel_coordinates
3605	//
3606	// @desc generate the pixel center coordinates
3607	//
3608	// @var *pcam structure camera containing all the
3609	// camera information
3610	// @return total number of pixels
3611	//
3612	// DP
3613	//
3614	// @date Thu Oct 14 10:41:03 CEST 1999
3615	//------------------------------------------------------------
3616	// @function
3617
3618	//!@{
3619	/****** igen_pixel_coordinates() *******************************/
3620
3621	int igen_pixel_coordinates(struct camera *pcam) {
3622	/* generate pixel coordinates, return value is number of pixels */
3623
3624	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
3625	float fsegment_fract;
3626	double dtsize;
3627	double dhsize;
3628	double dpsize;
3629	double dxfirst_pix;
3630	double dyfirst_pix;
3631	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5;
3632	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5;
3633
3634
3635	double dstartx, dstarty; /* for the gap pixels and outer pixels */
3636	int j, nrow;
3637
3638	dpsize = pcam->dpixdiameter_cm;
3639	dtsize = dpsize * sqrt(3.) / 2.;
3640	dhsize = dpsize / 2.;
3641
3642	/* Loop over central pixels to generate co-ordinates */
3643
3644	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
3645
3646	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
3647
3648	pcam->dpixsizefactor[ipixno-1] = 1.;
3649
3650	in = 0;
3651
3652	i = 0;
3653	itot_inside_ring = 0;
3654	iring_no = 0;
3655
3656	/* Calculate the number of pixels out to and including the ring containing pixel number */
3657	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
3658
3659	while (itot_inside_ring == 0){
3660
3661	iN = 3(i(i+1)) + 1;
3662
3663	if (ipixno <= iN){
3664	iring_no = i;
3665	itot_inside_ring = iN;
3666	}
3667
3668	i++;
3669	}
3670
3671
3672	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
3673
3674	ipix_in_ring = 0;
3675	for (i = 0; i < iring_no; ++i){
3676
3677	ipix_in_ring = ipix_in_ring + 6;
3678	}
3679
3680	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
3681	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
3682	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
3683
3684	isegment = 0;
3685	fsegment_fract = 0.;
3686	if (iring_no > 0) {
3687
3688	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
3689
3690	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
3691
3692	}
3693
3694	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
3695	/* pixel width (flat to flat) dpsize. */
3696
3697	dxfirst_pix = dpsize*iring_no;
3698	dyfirst_pix = 0.;
3699
3700	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
3701
3702	ddxseg1 = - dhsize*iring_no;
3703	ddyseg1 = dtsize*iring_no;
3704	ddxseg2 = -dpsize*iring_no;
3705	ddyseg2 = 0.;
3706	ddxseg3 = ddxseg1;
3707	ddyseg3 = -ddyseg1;
3708	ddxseg4 = -ddxseg1;
3709	ddyseg4 = -ddyseg1;
3710	ddxseg5 = -ddxseg2;
3711	ddyseg5 = 0.;
3712
3713	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
3714	/* anti-clockwise around the ring by adding the segment to segment vectors. */
3715
3716	switch (isegment) {
3717
3718	case 0:
3719
3720	pcam->dxc[ipixno-1] = 0.;
3721	pcam->dyc[ipixno-1] = 0.;
3722
3723	case 1:
3724	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
3725	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
3726
3727	break;
3728
3729	case 2:
3730
3731	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
3732	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
3733
3734	break;
3735
3736	case 3:
3737
3738	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
3739	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
3740
3741	break;
3742
3743	case 4:
3744
3745	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
3746	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
3747
3748	break;
3749
3750	case 5:
3751
3752	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
3753	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
3754
3755	break;
3756
3757	case 6:
3758
3759	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
3760	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
3761
3762	break;
3763
3764	default:
3765
3766	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
3767	return(0);
3768
3769	} /* end switch */
3770
3771	} /* end for */
3772
3773	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
3774	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
3775
3776	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
3777
3778	j = pcam->inumcentralpixels;
3779
3780	for(i=0; i<pcam->inumgappixels; i=i+6){
3781	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
3782	pcam->dyc[j + i ] = dstarty;
3783	pcam->dpixsizefactor[j + i] = 1.;
3784	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
3785	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
3786	pcam->dpixsizefactor[j + i + 1] = 1.;
3787	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
3788	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
3789	pcam->dpixsizefactor[j + i+ 2] = 1.;
3790	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
3791	pcam->dyc[j + i + 3] = dstarty;
3792	pcam->dpixsizefactor[j + i+ 3] = 1.;
3793	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
3794	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
3795	pcam->dpixsizefactor[j + i+ 4] = 1.;
3796	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
3797	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
3798	pcam->dpixsizefactor[j + i + 5] = 1.;
3799	} /* end for */
3800	} /* end if */
3801
3802	/* generate positions of the outer pixels */
3803
3804	if( pcam->inumbigpixels > 0 ){
3805
3806	j = pcam->inumcentralpixels + pcam->inumgappixels;
3807
3808	for(i=0; i<pcam->inumbigpixels; i++){
3809	pcam->dpixsizefactor[j + i] = 2.;
3810	}
3811
3812	in = 0;
3813
3814	nrow = (int) ceil(dstartx / 2. / dpsize);
3815
3816	while(in < pcam->inumbigpixels){
3817
3818	pcam->dxc[j + in] = dstartx + dpsize;
3819	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
3820	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
3821	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
3822	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
3823	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
3824	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
3825	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
3826	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
3827	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
3828	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
3829	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
3830	for(i=1; i<nrow; i++){
3831	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
3832	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
3833	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
3834	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
3835	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
3836	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
3837	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
3838	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
3839	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
3840	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
3841	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
3842	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
3843	}
3844	in = in + 6 * nrow;
3845	dstartx = dstartx + 2. * dpsize;
3846	nrow = nrow + 1;
3847	} /* end while */
3848
3849	} /* end if */
3850
3851	return(pcam->inumpixels);
3852
3853	}
3854	//!@}
3855
3856	//!-----------------------------------------------------------
3857	// @name bpoint_is_in_pix
3858	//
3859	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
3860	//
3861	// @var *pcam structure camera containing all the
3862	// camera information
3863	// @var dx, dy point coordinates in centimeters
3864	// @var ipixnum pixel number (starting at 0)
3865	// @return TRUE if the point is inside the pixel, FALSE otherwise
3866	//
3867	// DP
3868	//
3869	// @date Thu Oct 14 16:59:04 CEST 1999
3870	//------------------------------------------------------------
3871	// @function
3872
3873	//!@{
3874
3875	/****** bpoint_is_in_pix() *************************************/
3876
3877	#define sqrt13 0.577350269 // = 1./sqrt(3.)
3878	#define sqrt3 1.732050807 // = sqrt(3.)
3879
3880	int bpoint_is_in_pix(double dx, double dy, MGeomCam *pgeomcam)
3881	{
3882	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
3883	/* the pixel is assumed to be a "closed set" */
3884
3885	/*
3886	a = length of one of the edges of one pixel,
3887	b = half the width of one pixel
3888	*/
3889
3890	const int numN = pgeomcam->GetNumPixels();
3891
3892	for (int i=0; i<numN; i++)
3893	{
3894	MGeomPix &pixel = (*pgeomcam)[i];
3895	const double b = pixel.GetD()/2;
3896	const double a = pixel.GetD()/sqrt3;
3897
3898	const double xx = dx - pixel.GetX();
3899	const double yy = dy - pixel.GetY();
3900
3901	if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) \|\|
3902	((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){
3903
3904	return i; // inside i
3905	}
3906
3907	// return -1; // outside
3908	}
3909
3910	return -1; // outside
3911	}
3912
3913	//!@}
3914
3915	//------------------------------------------------------------
3916	// @name dist_r_P
3917	//
3918	// @desc distance straight line r - point P
3919	//
3920	// @date Sat Jun 27 05:58:56 MET DST 1998
3921	// @function @code
3922	//------------------------------------------------------------
3923	// dist_r_P
3924	//
3925	// distance straight line r (x+tl, y+tm, z+t*n) to point P(a,b,c)
3926	//
3927	// We assume that vector (l, m, n) is normalized l^2+m^2+n^2 = 1
3928	//------------------------------------------------------------
3929
3930	float
3931	dist_r_P(float a, float b, float c,
3932	float l, float m, float n,
3933	float x, float y, float z)
3934	{
3935	return (
3936	sqrt((SQR((a-x)m-(b-y)l) +
3937	SQR((b-y)n-(c-z)m) +
3938	SQR((c-z)l-(a-x)n)))
3939	);
3940	}
3941
3942	//------------------------------------------------------------
3943	// @name check_reflector_file
3944	//
3945	// @desc check if a given reflector file has the right signature
3946	// @desc return TRUE or FALSE
3947	//
3948	// @date Mon Feb 14 16:44:21 CET 2000
3949	// @function @code
3950	//------------------------------------------------------------
3951
3952	int check_reflector_file(FILE *infile){
3953
3954	char sign[20]; // auxiliary variable
3955
3956	strcpy(sign, REFL_SIGNATURE_B);
3957
3958	fread( (char *)sign, strlen(REFL_SIGNATURE_B), 1, infile);
3959	if (strcmp(sign, REFL_SIGNATURE_A) == 0){
3960	fread( (char *)sign, 1, 1, infile);
3961	return 4;
3962	}
3963	else if (strcmp(sign, REFL_SIGNATURE_B) == 0){
3964	fread( (char *)sign, 1, 1, infile);
3965	return 5;
3966	}
3967	else if (strcmp(sign, REFL_SIGNATURE_C) == 0){
3968	// An empty bin has been removed and therefore we do not need to rezd it.
3969	return 6;
3970	}
3971	else {
3972	cout << "ERROR: Signature of .rfl file is not correct\n";
3973	cout << '"' << sign << '"' << '\n';
3974	cout << "should be: " << REFL_SIGNATURE_A <<" or "<< REFL_SIGNATURE_B <<" or " << REFL_SIGNATURE_C <<" or "<< '\n';
3975	return(FALSE);
3976	}
3977
3978
3979	}
3980
3981	//------------------------------------------------------------
3982	// @name lin_interpol
3983	//
3984	// @desc interpolate linearly between two points returning the
3985	// @desc y-value of the result
3986	//
3987	// @date Thu Feb 17 11:31:32 CET 2000
3988	// @function @code
3989	//------------------------------------------------------------
3990
3991	float lin_interpol( float x1, float y1, float x2, float y2, float x){
3992
3993	if( (x2 - x1)<=0. ){ // avoid division by zero, return average
3994	cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
3995	return((y1+y2)/2.);
3996	}
3997	else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
3998	return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
3999	}
4000	}
4001
4002	//------------------------------------------------------------
4003	// @name size_rotated
4004	//
4005	// @desc It rotates the NSB
4006	//
4007	// @date Tue Jan 7 17:09:25 CET 2003
4008	// @function @code
4009	//------------------------------------------------------------
4010
4011	int size_rotated(
4012	float *rotated,
4013	float *nsb,
4014	float rho)
4015	{
4016	int r=0;
4017	float size_rotated;
4018	float Num_Pixels;
4019	float PixNum=0;
4020	float rho_pixel;
4021	int j=0;
4022	int k=0;
4023
4024	for(int i=1; i<iMAXNUMPIX;i++){
4025	// Substituir per Int_t radius[iMAXNUMPIX]={0,1,1,1,1,1,1,2,2,2,2, ...}
4026	Num_Pixels=0;
4027	for (int ii=1; ii<17;ii++){
4028	if (Num_Pixels >= i){
4029	r=ii-1;
4030	break;
4031	}
4032
4033
4034
4035	if (ii<12){
4036	Num_Pixels+=ii*6;
4037	PixNum=6*ii;
4038	// size_rotated = (360/PixNum);
4039	//rho_pixel = rho/size_rotated;
4040	}
4041
4042
4043	else
4044	{
4045	Num_Pixels+=(ii-6)*6;
4046	PixNum=(ii-6)*6;
4047	// size_rotated = (360/PixNum);
4048	// rho_pixel = rho/size_rotated;
4049	}
4050	}
4051
4052
4053
4054	size_rotated = (360/PixNum);
4055	rho_pixel = rho/size_rotated;
4056
4057	//Buscar j i k que no canvin de r!!!
4058
4059	j=i-int(rho_pixel);
4060
4061	//cout<<"j inicial "<<j<<endl;
4062	int MINPIXinner[13]= {0,1,7,19,37,61,91,127,169,217,271,331,397};
4063	int MINPIXouter[5]={397,433,475,523,578};
4064
4065	if( r<12 ){
4066	if (j < MINPIXinner[r])
4067	{
4068	j = i+6 * r - int(rho_pixel);
4069	}
4070
4071	k=j-1;
4072	if (k < MINPIXinner[r])
4073	{
4074	k = MINPIXinner[r+1]-1;
4075	}
4076	}
4077	if( r > 11 && r < 16){
4078	if (j < MINPIXouter[r-12])
4079	{
4080	j = i + (r - 6) * 6 - int(rho_pixel);
4081
4082	}
4083
4084	k=j-1;
4085	if ( k < MINPIXouter[r-12])
4086	{
4087	k = MINPIXouter[r-11] - 1;
4088	}
4089	}
4090	rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k];
4091	}
4092	return (int(rho_pixel));
4093
4094	}
4095
4096
4097	//------------------------------------------------------------
4098	// @name produce_phes
4099	//
4100	// @desc read the photons of an event, pixelize them and simulate QE
4101	// @desc return various statistics and the array of Photoelectrons
4102	//
4103	// @date Mon Feb 14 16:44:21 CET 2000
4104	// @function @code
4105	//------------------------------------------------------------
4106
4107	int produce_phes( FILE *sp, // the input file
4108	class MGeomCam *camgeom, // the camera layout
4109	float minwl_nm, // the minimum accepted wavelength
4110	float maxwl_nm, // the maximum accepted wavelength
4111	class MTrigger *trigger, // the generated phes
4112	class MFadc *fadc,
4113	int *itotnphe, // total number of produced photoelectrons
4114	float *nphe, // number of photoelectrons in each pixel
4115	int *incph, // total number of cph within camera
4116	float *tmin_ns, // minimum arrival time of all phes
4117	float *tmax_ns, // maximum arrival time of all phes
4118	int ict, // Telescope that is being analised to get the right QE.
4119	float mirror_fraction, // Fraction of total mirror really present
4120	float fadc_jitter // Random shift (max 1 slice) of pulses in
4121	// FADC, to simulate FADC clock noise.
4122	){
4123
4124	static uint i;
4125	static int k, ipixnum;
4126	static float cx, cy, wl, qe;
4127	static float cw;
4128	static MCCphoton photon;
4129	static float **qept;
4130	static char flag[SIZE_OF_FLAGS + 1];
4131	static float radius_mm;
4132	static UInt_t seed = (UInt_t)(get_seeds(0)*get_seeds(1));
4133
4134	float time;
4135	float pmtjit;
4136
4137	// reset variables
4138
4139	for ( i=0; i<camgeom->GetNumPixels(); ++i )
4140	nphe[i] = 0.0;
4141
4142
4143	*itotnphe = 0;
4144	*incph = 0;
4145
4146	radius_mm = camgeom->GetMaxRadius();
4147
4148	TRandom random;
4149	random.SetSeed(seed);
4150
4151	float C1, C2, C3, rho;
4152
4153	//- - - - - - - - - - - - - - - - - - - - - - - - -
4154	// read photons and "map" them into the pixels
4155	//--------------------------------------------------
4156
4157	// initialize CPhoton
4158
4159	photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
4160
4161	// read the photons data
4162
4163	// loop over the photons
4164
4165	while (1) {
4166
4167	photon.read(sp);
4168	photon.get_flag(flag);
4169
4170	if (isA( flag, FLAG_END_OF_EVENT ))
4171	{
4172	fseek (sp, SIZE_OF_FLAGS-photon.mysize(), SEEK_CUR);
4173	break;
4174	}
4175
4176	// Check if photon is inside trigger time range
4177
4178	time = photon.get_t() ;
4179
4180	if (time-*tmin_ns>TOTAL_TRIGGER_TIME) {
4181	continue;
4182	}
4183
4184	//
4185	// Account for possibly missing mirrors, or lower reflectivity...
4186	// mirror_fraction is the fraction of the total mirror really
4187	// working.
4188	//
4189	if (mirror_fraction < 1.)
4190	if (RandomNumber > mirror_fraction)
4191	continue;
4192
4193	//
4194	// Pixelization
4195	//
4196
4197	cx = photon.get_x();
4198	cy = photon.get_y();
4199
4200	if (Spotsigma > 0.)
4201	{
4202	cx += random.Gaus(0.,Spot_x);
4203	cy += random.Gaus(0.,Spot_y);
4204	}
4205	if (missPointing > 0.)
4206	{
4207	// We should take intot acount the rotation of the FoV
4208	C1 = 0.48 * sin(Zenith) - 0.87 * cos(Zenith) * cos(Azimutal);
4209	C3 = (0.87 * cos(Zenith) - 0.48 * sin(Zenith) * cos(Azimutal));
4210	C2 = sqrt( sin(Zenith) * sin(Zenith) * sin(Azimutal) * sin(Azimutal) + C3 * C3 );
4211	rho = acos( C1/C2 );
4212	rho=(sin(Azimutal)<0 ? (2 * 3.14159 - rho) : rho);
4213
4214	rho = rho*180/3.14159;
4215
4216	cx += (missP_xcos(rho)-missP_ysin(rho))/(10*camgeom->GetConvMm2Deg());
4217	cy += (missP_xsin(rho)+missP_ycos(rho))/(10*camgeom->GetConvMm2Deg());
4218	}
4219
4220
4221	// get wavelength
4222
4223	wl = photon.get_wl();
4224
4225	// cout << "wl " << wl << " radius " << sqrt(cxcx + cycy) << "\n";
4226
4227	// check if photon has valid wavelength and is inside outermost camera radius
4228
4229	if( (wl > maxwl_nm) \|\| (wl < minwl_nm) \|\| (sqrt(cxcx + cycy)*10 > radius_mm ) ){
4230	continue;
4231
4232	}
4233
4234	ipixnum = bpoint_is_in_pix(cx10, cy10, camgeom);
4235
4236	// -1 = the photon is in none of the pixels
4237	// 0 = the phton is in the central pixel, which is not used for trigger
4238	if (ipixnum==-1 \|\| ipixnum==0) {
4239	continue;
4240	}
4241
4242	// AM changed meaning of incph: before it was all photons read from
4243	// reflector file, now only those within a valid pixel:
4244	//
4245	// increase number of photons
4246	(*incph)++;
4247
4248	//+++
4249	// QE simulation
4250	//---
4251
4252	// set pointer to the QE table of the relevant pixel
4253
4254	qept = (float **)QE[ict][ipixnum];
4255
4256	// check if wl is inside table; outside the table, QE is assumed to be zero
4257
4258	if((wl < qept[0][0]) \|\| (wl > qept[0][pointsQE[ict]-1])){
4259	continue;
4260
4261	}
4262
4263	// find data point in the QE table (-> k)
4264
4265	k = 1; // start at 1 because the condition was already tested for 0
4266	while (k < pointsQE[ict]-1 && qept[0][k] < wl){
4267	k++;
4268	}
4269
4270	// calculate the qe value between 0. and 1.
4271
4272	qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
4273
4274	//
4275	// Apply incident angular correction due to Light Guides, plexiglas,
4276	// 1st dynode collection efficiency, double crossings... etc.
4277	// This information is contained in the file Data/LightCollection.dat
4278	//
4279	cw=photon.get_phi()*180./3.14159265;
4280
4281	// find data point in the WC table (-> k)
4282
4283	if ( camgeom->GetPixRatio(ipixnum) < 1.) // => Pixel is outer pixel
4284	{
4285	k = 0;
4286	while (k < pointsWC-1 && WC_outer[0][k] < cw)
4287	k++;
4288	// correct the qe with WC data.
4289	qe = qe*lin_interpol(WC_outer[0][k-1], WC_outer[1][k-1],
4290	WC_outer[0][k], WC_outer[1][k], cw);
4291	}
4292
4293	else // => Pixel is inner pixel
4294	{
4295	k = 0;
4296	while (k < pointsWC-1 && WC[0][k] < cw)
4297	k++;
4298	// correct the qe with WC data.
4299	qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
4300	}
4301
4302
4303	// if random > quantum efficiency, reject it
4304
4305	if ( (RandomNumber) > qe ) {
4306	continue;
4307	}
4308
4309	//+++
4310	// The photon has produced a photo electron
4311	//---
4312
4313	// cout << " accepted\n";
4314
4315	// increment the number of photoelectrons in the relevant pixel
4316
4317	nphe[ipixnum] += 1.0;
4318
4319	//
4320	// PMT time jitter: gaussian, not negative (MTrigger::FillShow does
4321	// not accept negative times!)
4322	//
4323	do
4324	pmtjit = random.Gaus(3.*pmt_jitter, pmt_jitter);
4325	while(pmtjit<0.);
4326
4327	// store the new photoelectron
4328
4329	fadc->Fill(ipixnum,
4330	(time-*tmin_ns) + pmtjit + fadc_jitter,
4331	trigger->FillShow(ipixnum, time-*tmin_ns + pmtjit + fadc_jitter),
4332	!((camgeom)[ipixnum].GetD()>(camgeom)[0].GetD()));
4333
4334	*itotnphe += 1;
4335	}
4336
4337	seed = random.GetSeed(); // Get seed for next call
4338
4339	return(0);
4340
4341	}
4342
4343
4344	//------------------------------------------------------------
4345	// @name produce_nsbrates
4346	//
4347	// @desc read the starfield file, call produce_phes on it in,
4348	// @desc different wavebands, calculate the nsbrates
4349	//
4350	// @date Mon Feb 14 16:44:21 CET 2000
4351	// @function @code
4352	//------------------------------------------------------------
4353
4354	int produce_nsbrates( char *iname, // the starfield input file name
4355	MGeomCam *camgeom, // camera layout
4356	float **rate_phepns,
4357	// the product of this function:
4358	// the NSB rates in phe/ns for each pixel
4359	int ict,
4360	float mirror_fraction)
4361	{
4362	uint i, j; // counters
4363	int k, ii; // counters
4364
4365	static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
4366	WAVEBANDBOUND2,
4367	WAVEBANDBOUND3,
4368	WAVEBANDBOUND4,
4369	WAVEBANDBOUND5,
4370	WAVEBANDBOUND6 };
4371
4372	FILE *infile; // the input file
4373	fpos_t fileposition; // marker on the input file
4374	static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
4375	static MCEventHeader evth; // the event header
4376	static MCEventHeader evth_2; // the event header
4377	static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
4378	int reflector_file_version;
4379	int itnphe; // total number of produced photoelectrons
4380	int itotnphe; // total number of produced photoelectrons after averaging
4381	int incph; // total number of cph read
4382	float tmin_ns; // minimum arrival time of all phes
4383	float tmax_ns; // maximum arrival time of all phes
4384	float integtime_ns; // integration time from the starfield generator
4385	char flag_new[4];
4386
4387
4388	if (strlen(iname) == 0)
4389	{
4390	log( SIGNATURE, "No starfield input file has been provided.\n");
4391	return (0);
4392	}
4393	else // check if the starfield input file exists and open it
4394	{
4395	log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
4396	infile = fopen( iname, "r" );
4397
4398	if ( infile == NULL )
4399	{
4400	log( SIGNATURE, "ERROR! Cannot open starfield input file: %s\n", iname );
4401	exit(-1);
4402	}
4403	}
4404
4405
4406	// get signature, and check it
4407
4408	if((reflector_file_version=check_reflector_file( infile ))==FALSE){
4409	exit(1);
4410	}
4411
4412	// Instance of MTrigger and MFadc; needed here only as dummies for
4413	// a call to produce_phes (see below).
4414	// 15/09/2004, A. Moralejo, changed "trigger" and "flashadc" to
4415	// pointers, the former static allocation caused memory problems and
4416	// segmentation fault in some systems.
4417
4418	MTrigger* trigger = new MTrigger(camgeom->GetNumPixels(),camgeom,
4419	Trigger_gate_length,
4420	Trigger_overlaping_time,
4421	Trigger_response_ampl,
4422	Trigger_response_fwhm);
4423
4424	MFadc* flashadc = new MFadc(camgeom->GetNumPixels(),
4425	FADC_shape,
4426	FADC_response_integ,FADC_response_fwhm,
4427	FADC_shape_out,
4428	FADC_resp_integ_out,FADC_resp_fwhm_out,
4429	get_trig_delay(),
4430	FADC_slices_per_ns,
4431	FADC_slices_written);
4432
4433	// initialize flag
4434	strcpy( cflag, " \0" );
4435
4436	// get flag
4437
4438	fread( cflag, SIZE_OF_FLAGS, 1, infile );
4439
4440	if ( ! feof(infile)){
4441
4442	// reading .rfl file
4443
4444	if(!isA( cflag, FLAG_START_OF_RUN )){
4445	error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
4446	}
4447	else { // found start of run
4448
4449	if (reflector_file_version > 5){
4450
4451	fread( flag_new, 4, 1, infile );
4452
4453	if(!isA( flag_new, FLAG_START_OF_HEADER)){
4454
4455	// We break the main loop
4456	cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
4457	cout<<" We assume no Star Light"<<endl;
4458	return(0);
4459	}
4460
4461	Float_t flag_temp[SIZE_OF_MCRUNHEADER];
4462
4463	fread( flag_temp, (SIZE_OF_MCRUNHEADER)*sizeof(float), 1, infile );
4464
4465	}
4466
4467	fread( cflag, SIZE_OF_FLAGS, 1, infile );
4468
4469	if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
4470	fread( flag_new, 4, 1, infile );
4471
4472	if(!isA( flag_new, FLAG_EVENT_HEADER)){
4473
4474	// We break while events loop
4475	cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<endl;
4476	cout<<" We assume no light from Stars"<<endl;
4477	return(0);
4478	}
4479
4480	// get MCEventHeader
4481
4482	if (reflector_file_version<6)
4483	// fread( (char*)&evth, evth.mysize(), 1, infile );
4484	evth.read(infile);
4485	else
4486	// fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
4487	evth_2.read(infile);
4488
4489	if (reflector_file_version<6)
4490	integtime_ns = evth.get_energy();
4491	else
4492	integtime_ns = evth_2.get_energy();
4493
4494	// memorize where we are in the file
4495
4496	if (fgetpos( infile, &fileposition ) != 0){
4497	error( SIGNATURE, "Cannot position in file ...\n");
4498	}
4499
4500	// loop over the wavebands
4501
4502	for(i=0; i<iNUMWAVEBANDS; i++){
4503
4504	// initialize the rate array
4505
4506	for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4507	rate_phepns[j][i] = 0.;
4508	}
4509
4510	itotnphe = 0;
4511
4512	// read the photons and produce the photoelectrons
4513	// - in order to average over the QE simulation, call the
4514	// production function iNUMNSBPRODCALLS times
4515
4516	for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
4517
4518	// position the file pointer to the beginning of the photons
4519
4520	fsetpos( infile, &fileposition);
4521
4522	// produce photoelectrons (photons from starfield)
4523
4524	k = produce_phes( infile,
4525	camgeom,
4526	wl_nm[i],
4527	wl_nm[i+1],
4528	trigger, // this is a dummy here
4529	flashadc, // this is a dummy here
4530	&itnphe,
4531	nphe, // we want this!
4532	&incph,
4533	&tmin_ns,
4534	&tmax_ns,
4535	0,
4536	mirror_fraction,
4537	0.);
4538
4539
4540	if( k != 0 ){ // non-zero returnvalue means error
4541	cout << "Exiting.\n";
4542	exit(1);
4543	}
4544
4545	for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4546	rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
4547	}
4548
4549	itotnphe += itnphe;
4550
4551	} // end for(ii=0 ...
4552
4553	fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
4554	wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
4555
4556	} // end for(i=0 ...
4557
4558	}
4559	else{
4560	cout << "Starfield file contains no event.\nExiting.\n";
4561	exit(1);
4562	} // end if( isA ... event
4563	} // end if ( isA ... run
4564	}
4565	else{
4566	cout << "Starfield file contains no run.\nExiting.\n";
4567	exit(1);
4568	}
4569
4570	fclose( infile );
4571
4572	return(0);
4573	}
4574
4575
4576
4577	//-------------------------------------------------------------
4578	//
4579	// Function DoCalibration. A. Moralejo October 2004
4580	//
4581	// Generates calibration events similar to those in the real
4582	// calibration runs.
4583	//
4584	//-------------------------------------------------------------
4585
4586	int DoCalibration(MFadc Fadc_CT, MTrigger Trigger_CT,
4587	TObjArray camgeom,
4588	float nsb_trigresp, float nsb_fadcresp,
4589	MLons lons, MLons lons_outer,
4590	float **nsb_phepns, int addElecNoise,
4591	TTree EvtTree, MRawEvtHeader *EvtHeader,
4592	MMcEvt McEvt, MRawEvtData EvtData)
4593	{
4594
4595	int nevent = 0;
4596	int maxevents;
4597	int lons_return;
4598
4599	int *ntotphe = new int[ct_Number];
4600	int *ntotphot = new int[ct_Number];
4601	float *nphe_from_nsb = new float[ct_Number];
4602	float *pulsetime = new float[ct_Number];
4603
4604	float lambda, sigma_lambda, phot_per_pix, sigma_time;
4605	UInt_t first_pixel, last_pixel;
4606	int selected_pixel;
4607
4608	TRandom random;
4609	random.SetSeed((UInt_t)(get_seeds(0)*get_seeds(0)));
4610
4611	for (int ict = 0; ict < ct_Number; ict++)
4612	pulsetime[ict] = 0.;
4613
4614	// Get parameters of the calibration:
4615
4616	get_calibration_properties( &lambda, &sigma_lambda, &phot_per_pix,
4617	&sigma_time, &maxevents, &selected_pixel);
4618
4619	first_pixel = new UInt_t[ct_Number];
4620	last_pixel = new UInt_t[ct_Number];
4621
4622	for (int ict = 0; ict < ct_Number; ict++)
4623	{
4624	if (selected_pixel < 0)
4625	{
4626	first_pixel[ict] = 0;
4627	last_pixel[ict] = ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels()-1;
4628	}
4629	else
4630	first_pixel[ict] = last_pixel[ict] = selected_pixel;
4631	}
4632
4633
4634	TArrayC *fadcValues; //@< the analog Fadc High gain signal for pixels
4635	TArrayC *fadcValuesLow; //@< the analog Fadc Low gain signal for pixels
4636
4637	// allocate space for FADC info
4638	fadcValues = new TArrayC(FADC_slices_written);
4639	fadcValuesLow = new TArrayC(FADC_slices_written);
4640
4641
4642	// allocate space for PMTs numbers of pixels
4643	float *pheinpix = new float [ct_NPixels];
4644
4645	for (int calevent = 0; calevent < maxevents; calevent++)
4646	{
4647	//
4648	// Clear Trigger and Fadc
4649	//
4650	for(int ict = 0; ict < ct_Number; ict++)
4651	{
4652	Trigger_CT[ict]->Reset() ;
4653	Trigger_CT[ict]->ClearFirst();
4654	Trigger_CT[ict]->ClearZero();
4655	Fadc_CT[ict]->Reset() ;
4656
4657	ntotphe[ict] = ntotphot[ict] = 0;
4658	nphe_from_nsb[ict] = 0.;
4659
4660	}
4661
4662	nevent++;
4663
4664	if((nevent+1)%100 == 1)
4665	log(SIGNATURE, "Event %d\n", nevent);
4666
4667	// Produce the photoelectrons
4668	for(int ict = 0; ict < ct_Number; ict++)
4669	{
4670	// Obtain the FADC jitter of 1 FADC slice. This is a time to be added to the
4671	// time of all photons in an event, before digitalization of the signal. It is
4672	// therefore the same time shift for all pixels in a CT.
4673
4674	float fadc_jitter =
4675	(1./Fadc_CT[ict]->GetFadcSlicesPerNanosec()) * random.Uniform(0., 1.); //ns
4676
4677
4678	produce_calib_phes( (MGeomCam*)(camgeom.UncheckedAt(ict)),
4679	Trigger_CT[ict],
4680	Fadc_CT[ict],
4681	&(ntotphe[ict]),
4682	pheinpix,
4683	&(ntotphot[ict]),
4684	ict,
4685	lambda,
4686	sigma_lambda,
4687	phot_per_pix,
4688	sigma_time,
4689	selected_pixel,
4690	fadc_jitter
4691	);
4692
4693	pulsetime[ict] = fadc_jitter; // Keep value for writing it to output.
4694	}
4695
4696	// NSB simulation
4697
4698	if(lons && lons_outer)
4699	{
4700	// Fill trigger and fadc response in the trigger class from the NSB database
4701
4702	for (int ict = 0; ict < ct_Number; ict++)
4703	{
4704	for (UInt_t ui = first_pixel[ict]; ui <= last_pixel[ict]; ui++)
4705	{
4706	nphe_from_nsb[ict] += nsb_phepns[ict][ui];
4707
4708	if (nsb_phepns[ict][ui] > 0.0)
4709	{
4710	if((((MGeomCam)(camgeom.UncheckedAt(ict))))[ui].GetD() >
4711	(((MGeomCam)(camgeom.UncheckedAt(ict))))[0].GetD())
4712	lons_return = lons_outer->GetResponse(nsb_phepns[ict][ui],0.01,
4713	& nsb_trigresp[0],
4714	& nsb_fadcresp[0]);
4715	else
4716	lons_return = lons->GetResponse(nsb_phepns[ict][ui],0.01,
4717	& nsb_trigresp[0],
4718	& nsb_fadcresp[0]);
4719
4720	if (lons_return == 0)
4721	{
4722	cout << "Exiting.\n";
4723	exit(1);
4724	}
4725
4726	Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
4727	}
4728	}
4729	}
4730
4731	}// end if(simulateNSB) ...
4732
4733
4734	//++++++++++++++++++++++++++++++++++++++++++++++++++
4735	// at this point we have a camera full of
4736	// ph.e.s
4737	//--------------------------------------------------
4738
4739	// now the noise of the electronic
4740	// (preamps, optical transmission,..) is introduced.
4741	//
4742
4743	for (int ict = 0; ict < ct_Number; ict++)
4744	{
4745	if (addElecNoise)
4746	Fadc_CT[ict]->ElecNoise();
4747
4748	// now a shift in the fadc signal due to the pedestals is
4749	// introduced
4750	// This is done inside the class MFadc by the method Pedestals
4751
4752	Fadc_CT[ict]->Pedestals();
4753
4754	// Set the trigger time. The 3 ns are such that the calibration pulses
4755	// appear roughly at the same position as for the case of real data.
4756	// If you want to shift the pulse position, do not change this value here,
4757	// use the option trigger_delay in the input card instead.
4758	// The additional value 3*sigma_time makes that the pulse maximum is, in
4759	// average, in the same position no matter of the time width of the pulse
4760	// (see that, in produce_calib_phes(...), in order to avoid negative times
4761	// we shift the time of the photons by this same amount!)
4762
4763	Fadc_CT[ict]->TriggeredFadc(3.+3*sigma_time);
4764
4765	// Add the "digital noise": electronic noise intrinsic to the FADC and which
4766	// therefore is not scaled down in the low gain slices!
4767
4768	Fadc_CT[ict]->DigitalNoise();
4769
4770	//
4771	// Fill the event header information
4772	//
4773	EvtHeader[ict]->FillHeader((UInt_t) calevent, 0);
4774
4775	// Set flag to indicate this is a calibration or a pedestal run:
4776	if (phot_per_pix > 1.e-3)
4777	{
4778	// Set trigger pattern. WARNING: the trigger pattern of real
4779	// data is inverted for the MRawRunHeader.fFormatVersion=5
4780	// which is the one we now set in the camera outputs. Hence
4781	// we use the binary inversion operators "~" below. The Calibration
4782	// pattern, on the contrary, is not inverted.
4783	EvtHeader[ict]->SetTriggerPattern((UInt_t)(~(MTriggerPattern::kCalibration \|
4784	MTriggerPattern::kTriggerLvl1)));
4785	//
4786	// FIXME! For now color and intensity of the pulser is fixed!
4787	EvtHeader[ict]->SetCalibrationPattern((UInt_t)((BIT(11) \| BIT(12))<<16));
4788	}
4789	else // 0 cal. photons per pixel ==> pedestal run
4790	EvtHeader[ict]->SetTriggerPattern((UInt_t)~MTriggerPattern::kPedestal);
4791
4792
4793	McEvt[ict]->Fill( 0, 0, 0., -1.0, -1.0, -1.0, 0., 0., 0., 0., 0.,
4794	0., 0., 0., 0., 0., 0., 0., 0.,
4795	0., 0., 0., 0., 0, 0, 0,
4796	(UInt_t) ntotphot[ict],
4797	(UInt_t) ntotphe[ict],
4798	(UInt_t) nphe_from_nsb[ict]+ ntotphe[ict],
4799	0., 0., 0., pulsetime[ict]);
4800
4801	//
4802	// Fill the FADC information
4803	//
4804	for(UInt_t ui = first_pixel[ict]; ui <= last_pixel[ict]; ui++)
4805	{
4806	for (int jslice = 0; jslice < FADC_slices_written; jslice++)
4807	{
4808	fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(ui, jslice),jslice);
4809	fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(ui,jslice),jslice);
4810	}
4811	EvtData[ict]->AddPixel(ui,fadcValues,0);
4812	EvtData[ict]->AddPixel(ui,fadcValuesLow,kTRUE);
4813	}
4814	}
4815
4816	EvtTree->Fill();
4817
4818	// clear all
4819	for(int ict=0;ict<ct_Number;ict++)
4820	{
4821	EvtHeader[ict]->Clear() ;
4822	EvtData[ict]->ResetPixels(0,0);
4823	McEvt[ict]->Clear() ;
4824	}
4825	}
4826
4827	return(0);
4828	}
4829
4830
4831	//------------------------------------------------------------
4832	//
4833	// Function produce_calib_phes, A. Moralejo Oct 2004
4834	//
4835	//------------------------------------------------------------
4836
4837	int produce_calib_phes( MGeomCam *camgeom, // The camera layout
4838	MTrigger *trigger,
4839	MFadc *fadc,
4840	int *itotnphe, // total number of produced photoelectrons
4841	float *nphe, // number of photoelectrons in each pixel
4842	int *nphot, // total number of photons in all pixels
4843	int ict, // Telescope that is being analised to get the right QE.
4844	float lambda, // Mean wavelength of light in nm
4845	float sigma_lambda, // Sigma of wavelengtgaussian spread
4846	float phot_per_pix, // Average # of photons per inner pixel
4847	float sigma_time, // Sigma of time spread of photons
4848	int selected_pixel,// if >= 0, only this pixel is used!
4849	float fadc_jitter // Random shift (max 1 slice) of pulses in
4850	// FADC, to simulate FADC clock noise.
4851	)
4852	{
4853
4854	int ipixnum;
4855	float cx, cy, wl, time, phi, phi_deg, qe;
4856	float **qept;
4857	float radius_mm, focal_dist_mm;
4858	int total_photons;
4859	float pmtjit;
4860
4861	static UInt_t seed = (UInt_t)(get_seeds(0)*get_seeds(1));
4862
4863	// reset variables
4864
4865	for ( uint i = 0; i < camgeom->GetNumPixels(); i++ )
4866	nphe[i] = 0.0;
4867
4868	*itotnphe = 0;
4869	*nphot = 0;
4870
4871	TRandom random;
4872	random.SetSeed(seed);
4873
4874	//
4875	// Create photons and "map" them into the pixels
4876	//
4877
4878	// Maximum radius of camera:
4879	radius_mm = camgeom->GetMaxRadius();
4880
4881	// Distance from center of mirror dish to camera plane:
4882	focal_dist_mm = camgeom->GetCameraDist()*1000;
4883
4884	// Cosine of the angle between telescope axis and line from center of mirror
4885	// dish to the edge of the camera:
4886	float cos_epsilon_max = cos(atan2(radius_mm, focal_dist_mm));
4887
4888
4889	// Calculate total number of photons to be produced.
4890	if (selected_pixel < 0)
4891	total_photons = (int) (phot_per_pix * 3.14159265 * radius_mm * radius_mm /
4892	(*camgeom)[0].GetA());
4893	else
4894	total_photons = (int) (phot_per_pix *
4895	(camgeom)[selected_pixel].GetA() / (camgeom)[0].GetA());
4896
4897	// loop over the photons
4898
4899	for (int iph = 0; iph < total_photons; iph++)
4900	{
4901	//
4902	// Simulate arrival times of photons to camera plane. We do not simulate the small
4903	// delays between pixels due to the different distances to the pulser.
4904	//
4905	// We do not want negative times, so center the gaussian at 3 sigma
4906	// and reject negative values:
4907
4908	do
4909	time = random.Gaus(3*sigma_time, sigma_time);
4910	while (time < 0.);
4911
4912	// wavelength
4913	wl = random.Gaus(lambda, sigma_lambda);
4914
4915	if( (wl > WAVEBANDBOUND6) \|\| (wl < WAVEBANDBOUND1))
4916	continue;
4917
4918	if (selected_pixel < 0)
4919	{
4920	// Obtain photon coordinates on the camera. We assume a point source of light placed
4921	// in the center of the mirror dish.
4922
4923	// polar angle
4924	float psi = RandomNumber * 2 * 3.14159265;
4925	// angle between the telescope axis and the photon trajectory.
4926	float epsilon = acos(1.-RandomNumber*(1.-cos_epsilon_max));
4927	float tanepsilon = tan(epsilon);
4928
4929	cx = focal_dist_mmtanepsiloncos(psi); // mm
4930	cy = focal_dist_mmtanepsilonsin(psi); // mm
4931
4932	if (sqrt(cxcx + cycy) > radius_mm )
4933	continue;
4934
4935	// Angle between photon trajectory and camera plane:
4936	phi = 3.14159265/2.-epsilon; // rad
4937
4938	//
4939	// Pixelization
4940	//
4941	ipixnum = bpoint_is_in_pix(cx, cy, camgeom);
4942
4943	// -1 = the photon is in none of the pixels
4944	// 0 = the phton is in the central pixel, which is not used
4945	if (ipixnum==-1 \|\| ipixnum==0)
4946	continue;
4947	}
4948	else
4949	{
4950	// Angle between photon trajectory and camera plane:
4951	phi = atan2( focal_dist_mm,
4952	sqrt( (camgeom)[selected_pixel].GetX()(*camgeom)[selected_pixel].GetX()+
4953	(camgeom)[selected_pixel].GetY()(*camgeom)[selected_pixel].GetY()));
4954
4955	ipixnum = selected_pixel;
4956	}
4957
4958
4959	// increase number of photons within pixels
4960	*nphot += 1;
4961
4962	//
4963	// QE simulation
4964	//
4965	// set pointer to the QE table of the relevant pixel
4966
4967	qept = (float **)QE[ict][ipixnum];
4968
4969	// check if wl is inside table; outside the table, QE is assumed to be zero
4970
4971	if((wl < qept[0][0]) \|\| (wl > qept[0][pointsQE[ict]-1]))
4972	continue;
4973
4974
4975	// find data point in the QE table (-> k)
4976	int k = 1; // start at 1 because the condition was already tested for 0
4977	while (k < pointsQE[ict]-1 && qept[0][k] < wl)
4978	k++;
4979
4980	// calculate the qe value between 0. and 1.
4981
4982	qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
4983
4984
4985	//
4986	// Apply incident angular correction due to Light Guides, plexiglas,
4987	// 1st dynode collection efficiency, double crossings... etc.
4988	// This information is contained in the file Data/LightCollection.dat,
4989	// and has been read into the array WC (which stands for "Winston Cones")
4990	//
4991	phi_deg = phi*180./3.14159265;
4992
4993	// find data point in the WC table (-> k)
4994
4995	if ( camgeom->GetPixRatio(ipixnum) < 1.) // => Pixel is outer pixel
4996	{
4997	k = 0;
4998	while (k < pointsWC-1 && WC_outer[0][k] < phi_deg)
4999	k++;
5000	// correct the qe with WC data.
5001	qe = qe*lin_interpol(WC_outer[0][k-1], WC_outer[1][k-1],
5002	WC_outer[0][k], WC_outer[1][k], phi_deg);
5003	}
5004
5005	else // => Pixel is inner pixel
5006	{
5007	k = 0;
5008	while (k < pointsWC-1 && WC[0][k] < phi_deg)
5009	k++;
5010	// correct the qe with WC data.
5011	qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], phi_deg);
5012	}
5013
5014
5015	// if random > quantum efficiency, reject it
5016
5017	if ( (RandomNumber) > qe ) {
5018	continue;
5019	}
5020
5021	//
5022	// The photon has produced a photo electron
5023	//
5024
5025	// increment the number of photoelectrons in the relevant pixel
5026
5027	nphe[ipixnum] += 1.;
5028
5029	//
5030	// PMT time jitter: gaussian, not negative (MTrigger::FillShow does
5031	// not accept negative times!)
5032	//
5033	do
5034	pmtjit = random.Gaus(3.*pmt_jitter, pmt_jitter);
5035	while(pmtjit<0.);
5036
5037	// store the new photoelectron
5038
5039	fadc->Fill(ipixnum, time + pmtjit + fadc_jitter,
5040	trigger->FillShow(ipixnum, time + pmtjit + fadc_jitter),
5041	!((camgeom)[ipixnum].GetD()>(camgeom)[0].GetD()));
5042
5043	*itotnphe += 1;
5044	}
5045
5046	seed = random.GetSeed(); // Get seed for next call
5047
5048	return(0);
5049
5050	}
5051
5052	// @endcode
5053
5054
5055	//=------------------------------------------------------------
5056	//!@subsection Log of this file.
5057
5058	//!@{
5059	//
5060	// $Log: not supported by cvs2svn $
5061	// Revision 1.89 2005/02/17 09:15:28 moralejo
5062	//
5063	// Set as default option that of writing all event headers to output file,
5064	// not only those of the triggered events. To disable it, set the input card
5065	// flag "no_write_all_event_headers".
5066	//
5067	// Changed such that output images for events below the minimum number of
5068	// photoelectrons nphe2NSB required to simulate the noise (NSB & electronic)
5069	// will be empty. This will avoid the problem of these events being processed,
5070	// without any noise, later in the chain. Although those images are not in the
5071	// output, one can still check in the headers (MMcTrig) how many such events
5072	// with less than nphe2NSB photoelectrons would have triggered.
5073	//
5074	// Revision 1.88 2005/02/11 20:00:01 moralejo
5075	//
5076	// Added to output container "MMcEvtBasic" with the most important MC
5077	// parameters, to be kept for all events (triggered or not) through
5078	// the whole analysis chain in order to allow the calculation of effective
5079	// areas.
5080	//
5081	// Updated version of MRawRunHeader from 4 to 5. This means the trigger and
5082	// calibration patterns are correctly set and can be decoded.
5083	//
5084	// Revision 1.87 2005/02/10 19:28:10 moralejo
5085	//
5086	// Substituted input card option "write_all_events" by
5087	// "write_all_event_headers". If set, this will make camera to write out
5088	// the event headers of all the processed events, no matter whether they
5089	// have triggered or not.
5090	//
5091	// Revision 1.86 2005/02/10 14:49:21 moralejo
5092	//
5093	// Added setting of trigger flags in MRawEvtHeader::fTriggerPattern:
5094	// Lvl1 flag for normal shower runs, Calibration flag for calibration
5095	// runs and Pedestal flag for "calibration runs" in which the number
5096	// of photons per pixel is set to 0.
5097	//
5098	// Revision 1.85 2005/02/10 12:00:32 moralejo
5099	//
5100	// Changed call to EvtHeader[ict]->FillHeader in line 2535: second argument
5101	// was 20, and I set it now to 0. No idea why that was 20: like that, it was
5102	// setting a CL trigger pattern on all events (???!!!)
5103	//
5104	// Revision 1.84 2004/12/15 01:56:39 moralejo
5105	//
5106	// Added input card option pmt_jitter_ns to simulate the time jitter of
5107	// PMTs. The input parameter is the sigma of a gaussian, which by
5108	// default is sigma=0.25 ns. This jitter is applied to each phe-
5109	// independently. We have not applied this to the NSB noise, since the
5110	// arrival time of NSB photons is random and nothing would change.
5111	//
5112	// Revision 1.83 2004/11/17 11:43:13 moralejo
5113	//
5114	// Made the necessary changes in starresponse to account for the new option
5115	// to switch off gain fluctuations for the NSB noise database generation.
5116	// The option in the Star Response card is "gain_fluctuations_off". The
5117	// version number of the NSB database has been updated to 1004 (see
5118	// MStarLight.hxx), now including information on whether the gain
5119	// fluctuations were on or off when the NSB database was generated.
5120	//
5121	// Revision 1.82 2004/11/17 11:34:49 moralejo
5122	//
5123	// Added input card command "noise_gain_fluctuations_off", to disable the
5124	// PMT gain fluctuations also for the NSB noise (just for tests). Added
5125	// a flag to MMcFadcHeader and MMcTrigHeader about this. Also copied the flag
5126	// fGainFluctuations of MMcFadcHeader to MMcTrigHeader, since the gain
5127	// fluctuations affect both the trigger and the signal in the FADC.
5128	//
5129	// Revision 1.81 2004/11/04 22:00:51 moralejo
5130	//
5131	// Removed unused variables fTelesTheta and fTelesPhi from MMcRunHeader. They
5132	// were not useful because telescope orientation may change from event to
5133	// event. Added fMirrorFraction to the MMcConfigRunHeader container in the
5134	// camera output.
5135	//
5136	// Revision 1.80 2004/10/26 19:21:05 moralejo
5137	//
5138	// Added fFadcTimeJitter to root output, in container MMcEvt. Added
5139	// also fGainFluctuations boolean flag to MMcFadcHeader, to keep track of
5140	// whether PMT gain fluctuations are simulated or not.
5141	//
5142	// Revision 1.79 2004/10/26 14:02:32 moralejo
5143	//
5144	// Added option to switch off gain fluctuations (gain_fluctuations_off)
5145	//
5146	// Revision 1.78 2004/10/21 17:44:07 moralejo
5147	//
5148	// Fixed error recently introduced in MLons::GetResponse. The NSB database
5149	// was not correclt ycopied to the FADC when the randomly selected time was
5150	// too close to the end of the database. This happened about 2% of times and
5151	// produced some repetitive noise peaks in the FADC.
5152	// Changed in camera.cxx the arguments of lons.SetSeed y lons_outer.SetSeed
5153	//
5154	// Revision 1.77 2004/10/19 10:35:05 moralejo
5155	// * empty log message *
5156	//
5157	// Revision 1.76 2004/10/14 16:56:43 moralejo
5158	//
5159	// - Added calibration_run option to produce calibration MC files.
5160	//
5161	// - Added jitter of pulse position +- 0.5 slices due to FADC clock noise.
5162	//
5163	// Revision 1.75 2004/10/14 12:55:02 moralejo
5164	// * empty log message *
5165	//
5166	// Revision 1.74 2004/10/13 17:05:05 moralejo
5167	// * empty log message *
5168	//
5169	// Revision 1.72 2004/10/12 13:39:34 moralejo
5170	//
5171	// Lots of changes intended to make it possible to select the FADC sampling
5172	// frequency from the input card, through the command fadc_GHz. The most
5173	// important ones are the following:
5174	//
5175	// - Removed FADC_SLICES de Mars/mmc/Mdefine.h
5176	// Already defined in MFadcDefine.h!
5177	//
5178	// - Replaced fixed numbers in array dimensions of starresponse.cxx
5179	//
5180	// - Added in MFadc.cxx and MFadc.hxx (Float_t) casts to initializations
5181	// of single phe response array
5182	//
5183	// - IMPORTANT: Fixed MFadc::GetResponse -> The returned single phe response
5184	// had only RESPONSE_SLICES (which is actually for the trigger branch),
5185	// whereas it should have RESPONSE_SLICES_MFADC. Fixed the same confusion
5186	// in other two points of the code (filling of the FADC for the signal),
5187	// in Fill() and FillOuter().
5188	//
5189	// - RESPONSE_SLICES_MFADC is now eliminated, since this quantity is now
5190	// decided by MFadc depending on other parameters.
5191	//
5192	// - Fixed problem in starresponse.cxx due to which the histograms fadcresp
5193	// and fadcbase in the root output were actually identical.
5194	//
5195	// - Changed starresponse.cxx such that above 1 phe/ns/pixel the precision
5196	// of the database is forced to be 0.1, and below 1, to 0.01 phe/ns/pixel
5197	//
5198	// - In Camera/creadparam.cxx: Now unkown tokens cause the program to stop,
5199	// instead of being simply skipped as it was until now.
5200	//
5201	// - Fixed error in MStarLight::StoreHisto. TH1::SetBinContent uses bin numbers
5202	// from 1 to nbins (in old code 0 to nbins-1 was assumed).
5203	//
5204	// - In MLons.cxx: use memcpy to copy pieces of the database into the FADC and
5205	// trigger branches, to (hopefully) speed up execution. For this I had to add
5206	// 2 getter functions in MStarLight.hxx
5207	//
5208	// - Everywhere: changed the shape parameter for trigger and FADC response to
5209	// be an Int. Changed version of NSB database from 1002 to 1003.
5210	//
5211	// - In MTrigger.cxx, changed all initializations of SlicesFirst and
5212	// SlicesSecond to 0 (instead of -50 as it was before). This controls the
5213	// time of trigger. If no trigger happened (like when making pedestal files)
5214	// the time was negative and the array index used to retrieve the noise from
5215	// the database was negative, resulting in "discontinuities" in the noise
5216	// (half-photoelectrons for instance).
5217	//
5218	// - In MStarLight changed fBinsTrig and fBinsFadc from Float_t to Int_t
5219	//
5220	// - Replace WIDTH_FADC_TIMESLICE by FADC_SLICES_PER_NSEC (which is the
5221	// inverse of the former)
5222	//
5223	// - Replaced SLICES_PER_NSEC by TRIG_SLICES_PER_NSEC
5224	//
5225	// - TRIGBINS eliminated. It depends on other two values
5226	// TRIGBINS = TIMERANGE*TRIG_SLICES_PER_NSEC
5227	//
5228	// - FADCBINS eliminated. It depends on other two values
5229	// FADCBINS = TIMERANGE*FADC_SLICES_PER_NSEC
5230	//
5231	// - MTriggerDefine.h Changed RESPONSE_SLICES to RESPONSE_SLICES_TRIG
5232	//
5233	// - Added to the MLons constructor an argument regarding the FADC sampling
5234	// frequency
5235	//
5236	// - MFadc: now the number of response slices for the FADC simulation is
5237	// decided by the program according to the other parameters.
5238	//
5239	// - MStarLight: Added arguments (fadc_slices_per_ns, response_slices_fadc)
5240	// to constructor.
5241	//
5242	// - creadparam.cxx, camera.cxx Changed default value of digital_noise to 0.
5243	//
5244	// Revision 1.71 2004/09/17 09:20:52 moralejo
5245	//
5246	// Updated some calls to current version of Mars:
5247	//
5248	// - EvtData[i]->InitRead(RunHeader) instead of EvtData[i]->Init(RunHeader);
5249	// - MRawRunHeader::kMagicNumber instead of just kMagicNumber
5250	// - EvtData[i]->ResetPixels (0, 0) instead of EvtData[i]->DeletePixels();
5251	//
5252	// Revision 1.70 2004/09/16 15:23:12 moralejo
5253	//
5254	// Changed "flashadc" and "trigger" in procedure produce_nsbrates from
5255	// objects to pointers (followed by dynamical allocation). This is only
5256	// to avoid memory problems (-> segmentation fault) in some systems.
5257	// Introduced missing initialization to 0 of *itotnphe in produce_phes.
5258	// Now the number of phes produced by stars shown on the screen make sense.
5259	//
5260	//
5261	// Revision 1.69 2004/03/30
5262	// Changed calculation of MMcFadcHeader.fPedesSigmaHigh and
5263	// MMcFadcHeader.fPedesSigmaLow to do as in real data (see comments in
5264	// code). Changed meaning of MMcFadcHeader.fAmplFadc and fAmplFadcOuter,
5265	// from amplitude to integral of single photoelectron pulse in FADC
5266	// counts. Added possibility to choose a realistic pulse shaped (as
5267	// measured using pulpo). Changed file Data/lightguides.dat by
5268	// Data/LightCollection.dat, intended to contain the information on
5269	// light collection efficiency regarding Winston cones, plexiglas, double
5270	// PMT crossing and colection efficiency of 1st dynode of PMT. Now the
5271	// information for inner and outer pixels can be different, since in the
5272	// LightCollection.dat file they are set independently.
5273	//
5274	// Revision 1.68 2004/02/06 17:39:24 blanch
5275	// Compiling with root 3.05 and updating MARS files.
5276	//
5277	// Revision 1.67 2004/01/30 09:51:18 blanch
5278	// [Changes mainly done by A. Moralejo]
5279	//
5280	// Several minnor changes have been done. For instance, some name of the
5281	// variables have been modified to a more self-explained name and
5282	// modifications while reading the asciis files at the end of the reflector file
5283	// has been introduced.
5284	//
5285	// The possibilty to enlarge the point spread function has been introduced
5286	// in order to be able to simualte the current data.
5287	//
5288	// All pixels are always written.
5289	//
5290	// Revision 1.65 2003/10/26 19:43:00 blanch
5291	// - The screen output information has been improved to prevent some
5292	// non-desired running conditions just looking at the output screen.
5293	// - One MMcEvt for each Telsecope is stored in the output file.
5294	// - 500 empty events are simualted to get a more precise estimation of the
5295	// pedestal Sigma for each pixel.
5296	//
5297	// Revision 1.64 2003/10/21 07:42:50 blanch
5298	// A factor 2.35 to transform the fwhm into the sigma of gaussian was missing
5299	// in the storing of FADC single hpe pulse determination.
5300	//
5301	// Revision 1.63 2003/10/17 19:38:31 blanch
5302	// Now the camera program will stop if a undefined Geometry is required.
5303	// The NSB is internally scaled for any camera geometry and qe.
5304	// The seeds in the input card are used to initilize the random numbers.
5305	// The Amplitud stored in the MMcFadcHeader is the amplitud of the sphe reponse.
5306	// The Pedestal rms is simulated in an artificial empty event.
5307	//
5308	// Revision 1.62 2003/09/26 11:25:07 blanch
5309	// Modification to be able to read MGeomCam branch for any Geometry.
5310	//
5311	// Revision 1.61 2003/09/25 17:09:20 blanch
5312	// Bug on the number of phe from diffuse NSB fixed.
5313	//
5314	// Revision 1.60 2003/09/23 16:50:55 blanch
5315	// WE do not read ct_file anymore since all Telescope information is
5316	// in the reflector or in MGeomCam.
5317	//
5318	// Revision 1.58 2003/09/15 09:59:53 blanch
5319	// The concept of the camera prgoram has not changed but this version has
5320	// quite a lot of changes to allow several Camera geometries as well as
5321	// multitelescope.
5322	//
5323	// It is suposed to be the first working code for camera 0.7.
5324	//
5325	// Revision 1.57 2003/07/17 18:02:46 blanch
5326	// Several new features introduced as well as fixed bugs
5327	//
5328	// - 1/100 events printed out
5329	// - Low gain implemented
5330	// - Different response for outer and inner pixels
5331	// - Some warnings removed
5332	// - pedestal and qe file from inpuit card
5333	// - Faster electronic simulation
5334	//
5335	// Revision 1.55 2003/02/12 12:22:10 blanch
5336	// * empty log message *
5337	//
5338	// Revision 1.54 2003/02/12 11:55:01 blanch
5339	// * empty log message *
5340	//
5341	// Revision 1.53 2003/01/23 18:35:21 blanch
5342	// * empty log message *
5343	//
5344	// Revision 1.52 2003/01/20 17:19:20 blanch
5345	// It fills the MMcCorsikaRun.
5346	//
5347	// Revision 1.51 2003/01/14 13:40:17 blanch
5348	// MRawRunHeader::fNumEvents has been filled with the number of events in
5349	// this file.
5350	// Problems in fImpact computation have been solved.
5351	// Option to set a dc value to rise the discriminator threshold has been added.
5352	//
5353	// Revision 1.50 2003/01/07 16:33:31 blanch
5354	// Star Field Rotation has been implemented by Raul Orduna. Now there is a
5355	// rotation for each shower. It is done by a non enter pixel rotation assuming
5356	// a circular symetry of the camera. It is not exact but it is accurate enough and
5357	// much faster.
5358	//
5359	// Revision 1.49 2002/12/13 10:04:07 blanch
5360	// * empty log message *
5361	//
5362	// Revision 1.48 2002/12/12 17:40:50 blanch
5363	// * empty log message *
5364	//
5365	// Revision 1.47 2002/12/10 17:19:31 blanch
5366	// * empty log message *
5367	//
5368	// Revision 1.46 2002/11/08 17:51:00 blanch
5369	// * empty log message *
5370	//
5371	// Revision 1.45 2002/10/29 17:15:27 blanch
5372	// Reading several reflector versions.
5373	//
5374	// Revision 1.44 2002/10/18 16:53:03 blanch
5375	// Modification to read several reflector files.
5376	//
5377	// Revision 1.43 2002/09/13 10:53:39 blanch
5378	// Minor change to remove some undisired comments.
5379	//
5380	// Revision 1.42 2002/09/09 16:00:49 blanch
5381	// Statement has been included to avoid writting to disk MParContainer and MArray.
5382	// It has also been added the effect of the WC, the actual values must be added,
5383	// once they are measured.
5384	//
5385	// Revision 1.41 2002/09/04 09:57:42 blanch
5386	// Modifications done to use MGeomCam from MARS.
5387	//
5388	// Revision 1.40 2002/07/16 16:15:22 blanch
5389	// A first implementation of the Star field rotation has been done.
5390	//
5391	// Revision 1.39 2002/04/27 10:48:39 blanch
5392	// Some unused varibles have been removed.
5393	//
5394	// Revision 1.38 2002/03/18 18:44:29 blanch
5395	// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
5396	//
5397	// Revision 1.37 2002/03/18 16:42:20 blanch
5398	// The data member fProductionSite of the MMcRunHeader has been set to 0,
5399	// which means that the prodution site is unknown.
5400	//
5401	// Revision 1.35 2002/03/15 15:15:52 blanch
5402	// Several modifications to simulate the actual trigger zone.
5403	//
5404	// Revision 1.34 2002/03/13 18:13:56 blanch
5405	// Some changes to fill correctly the new format of MMcRunHeader.
5406	//
5407	// Revision 1.33 2002/03/04 17:21:48 blanch
5408	// Small and not important changes.
5409	//
5410	// Revision 1.32 2002/02/28 15:04:52 blanch
5411	// A small back has been solved. Before, while not using the option
5412	// writte_all_images, not all triggered showers were stored. Now it is solved.
5413	// For that it is important that the less restrictive trigger option is
5414	// checked first.
5415	// A new facility has been introduced and now one can choose the step size in
5416	// trigger loop mode for the discriminator threshold.
5417	// The close-compact topology for less than 3 pixels does not make sense. Before
5418	// the program was ignoring that, now it switch to simple neighbour condition.
5419	//
5420	// Revision 1.31 2002/01/18 17:41:02 blanch
5421	// The option of adding noise to all pixels or to not adding the noise
5422	// has been added.
5423	// We removed the pixels larger than 577. When there were more than one
5424	// trigger in one shower, the pixel number was increasing. Now it is
5425	// flagged by the variable MMcTrig::fFirstLvlTrig.
5426	//
5427	// Revision 1.30 2001/11/27 09:49:54 blanch
5428	// Fixing bug which was treating wrongly the extension of star photons.
5429	//
5430	// Revision 1.29 2001/11/14 17:38:23 blanch
5431	// Sveral changes have been done:
5432	// - bpoint_in_in_pixel has been dodified to speed up the program
5433	// - Some minor changes have been done to remove warnings
5434	// - buffer size and split version of the Branches have been removed
5435	// - Some modifications were needed to adat the program to the new
5436	// MRawEvtData::DeletePixels
5437	//
5438	// Revision 1.28 2001/10/26 16:31:45 blanch
5439	// Removing several warnings.
5440	//
5441	// Revision 1.27 2001/09/05 10:04:33 blanch
5442	// * empty log message *
5443	//
5444	// Revision 1.26 2001/07/19 09:29:53 blanch
5445	// Different threshold for each pixel can be used.
5446	//
5447	// Revision 1.25 2001/05/08 08:07:54 blanch
5448	// New numbering for branches from different trigger conditions has been
5449	// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
5450	// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
5451	// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
5452	// information already filled.
5453	//
5454	// Revision 1.24 2001/04/06 16:48:09 magicsol
5455	// New camera version able to read the new format of the reflector output:
5456	// reflector 0.4
5457	//
5458	// Revision 1.23 2001/03/28 16:13:41 blanch
5459	// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
5460	// been solved.
5461	//
5462	// Revision 1.22 2001/03/20 18:52:43 blanch
5463	// * empty log message *
5464	//
5465	// Revision 1.21 2001/03/19 19:53:03 blanch
5466	// Some print outs have been removed.
5467	//
5468	// Revision 1.20 2001/03/19 19:30:06 magicsol
5469	// Minor changes have been done to improve the FADC pedestals treatment.
5470	//
5471	// Revision 1.19 2001/03/05 11:14:41 magicsol
5472	// I changed the position of readinf a parameter. It is a minnor change.
5473	//
5474	// Revision 1.18 2001/03/05 10:36:52 blanch
5475	// A branch with information about the FADC simulation (MMcFadcHeader) is writen
5476	// in the McHeaders tree of the aoutput root file.
5477	// The NSB contribution is only applied if the the number of phe form the shower
5478	// are above a given number.
5479	//
5480	// Revision 1.17 2001/02/23 11:05:57 magicsol
5481	// Small changes due to slightly different output format and the introduction of
5482	// pedesals for teh FADC.
5483	//
5484	// Revision 1.16 2001/01/18 18:44:40 magicsol
5485	// * empty log message *
5486	//
5487	// Revision 1.15 2001/01/17 09:32:27 harald
5488	// The changes are neccessary to have the same name for trees in MC and in
5489	// data. So now there should be no differences in MC data and real data from
5490	// FADC system.
5491	//
5492	// Revision 1.14 2001/01/15 12:33:34 magicsol
5493	// Some modifications have been done to use the new (Dec'2000) Raw data format.
5494	// There are also some minnors modifications to adapt the improvements in the
5495	// MTrigger class (overlaping time and trigger cells).
5496	//
5497	// Revision 1.13 2000/10/25 08:14:23 magicsol
5498	// The routine that produce poisson random numbers to decide how many phe
5499	// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
5500	// is used.
5501	//
5502	// Revision 1.10 2000/07/04 14:10:20 MagicSol
5503	// Some changes have been done in the root output file. The RawEvt tree is only
5504	// stored in the single trigger mode.
5505	// The trigger input parameters are also given by the general input card.
5506	// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
5507	// each one can be studied seperately.
5508	//
5509	// Revision 1.9 2000/06/13 13:25:24 blanch
5510	// The multiple arrays have been replaced, since they do not work
5511	// in alpha machines. Now we are using pointers and the command new
5512	// to allocate memory.
5513	//
5514	// Revision 1.8 2000/05/11 13:57:27 blanch
5515	// The option to loop over several trigger configurations has been included.
5516	// Some non-sense with triggertime range has been solved.
5517	// Montecarlo information and ADC counts are saved in a root file.
5518	// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
5519	//
5520	// Revision 1.7 2000/03/24 18:10:46 blanch
5521	// A first FADC simulation and a trigger simulation are already implemented.
5522	// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
5523	// 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.
5524	//
5525	// Revision 1.6 2000/03/20 18:35:11 blanch
5526	// 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.
5527	//
5528	// Revision 1.5 2000/02/18 17:40:35 petry
5529	// This version includes drastic changes compared to camera.cxx 1.4.
5530	// It is not yet finished and not immediately useful because the
5531	// trigger simulation is not yet re-implemented. I had to take it
5532	// out together with some other stuff in order to tidy the whole
5533	// program up. This is not meant as an insult to anyone. I needed
5534	// to do this in order to be able to work on it.
5535	//
5536	// This version has been put in the repository in order to be
5537	// able to share the further development with others.
5538	//
5539	// If you need something working, wait or take an earlier one.
5540	// See file README.
5541	//
5542	// Revision 1.4 2000/01/25 08:36:23 petry
5543	// The pixelization in previous versions was buggy.
5544	// This is the first version with a correct pixelization.
5545	//
5546	// Revision 1.3 2000/01/20 18:22:17 petry
5547	// Found little bug which makes camera crash if it finds a photon
5548	// of invalid wavelength. This bug is now fixed and the range
5549	// of valid wavelengths extended to 290 - 800 nm.
5550	// This is in preparation for the NSB simulation to come.
5551	// Dirk
5552	//
5553	// Revision 1.2 1999/11/19 08:40:42 harald
5554	// Now it is possible to compile the camera programm under osf1.
5555	//
5556	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
5557	// This the starting point for CVS controlled further developments of the
5558	// camera program. The program was originally written by Jose Carlos.
5559	// But here you can find a "rootified" version to the program. This means
5560	// that there is no hbook stuff in it now. Also the output of the
5561	// program changed to the MagicRawDataFormat.
5562	//
5563	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
5564	//
5565	// Revision 1.3 1999/10/22 15:01:28 petry
5566	// version sent to H.K. and N.M. on Fri Oct 22 1999
5567	//
5568	// Revision 1.2 1999/10/22 09:44:23 petry
5569	// first synthesized version which compiles and runs without crashing;
5570	//
5571	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
5572	// first synthesised version
5573	//
5574	// Revision 1.13 1999/03/15 14:59:05 gonzalez
5575	// camera-1_1
5576	//
5577	// Revision 1.12 1999/03/02 09:56:10 gonzalez
5578	// * empty log message *
5579	//
5580	//
5581	//!@}
5582
5583	//=EOF

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

Download in other formats: