Context Navigation

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

Visit:

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

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

Download in other formats: