Context Navigation

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

Visit:

Last change on this file since 6545 was 6386, checked in by moralejo, 20 years ago
Added to output container "MMcEvtBasic" with the most important MC parameters, to be kept for all events (triggered or not) through the whole analysis chain in order to allow the calculation of effective areas. Updated version of MRawRunHeader from 4 to 5. This means the trigger and calibration patterns are correctly set and can be decoded.
File size: 160.9 KB

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

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

Download in other formats: