Context Navigation

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

Visit:

Last change on this file since 6365 was 6365, checked in by moralejo, 20 years ago
Substituted input card option "write_all_events" by "write_all_event_headers". If set, this will make camera to write out the event headers of all the processed events, no matter whether they have triggered or not.
File size: 159.4 KB

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

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

Download in other formats: