Context Navigation

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

Visit:

Last change on this file since 5331 was 5320, checked in by moralejo, 21 years ago
Added fFadcTimeJitter to root output, in container MMcEvt. Added also fGainFluctuations boolean flag to MMcFadcHeader, to keep track of whether PMT gain fluctuations are simulated or not.
File size: 154.5 KB

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

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

Download in other formats: