Context Navigation

camera.cxx@ 6739

Visit:

Last change on this file since 6739 was 6739, checked in by moralejo, 20 years ago
Change in Mars/mmc/MFadcDefine.h: changed the value of SUBBINS from 5 to 500 to allow the proper simulation of the PMT time jitter. The analog signal is simulated in time steps equal to the duration of a FADC slice divided by SUBBINS. So now the binning is 100 finer than before. This will slow down the program somehow.
File size: 166.0 KB

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

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

Context Navigation

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

Download in other formats: