Context Navigation

camera.cxx@ 5440

Visit:

Last change on this file since 5440 was 5419, checked in by moralejo, 21 years ago

Made the necessary changes in starresponse to account for the new option to switch off gain fluctuations for the NSB noise database generation. The option in the Star Response card is "gain_fluctuations_off". The version number of the NSB database has been updated to 1004 (see MStarLight.hxx), now including information on whether the gain fluctuations were on or off when the NSB database was generated.

File size: 157.2 KB

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

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

Context Navigation

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

Download in other formats: