Context Navigation

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

Visit:

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

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

Download in other formats: