Context Navigation

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

Visit:

Last change on this file since 6351 was 6351, checked in by moralejo, 20 years ago
Added setting of trigger flags in MRawEvtHeader::fTriggerPattern: Lvl1 flag for normal shower events, Calibration flag for calibration events and Pedestal flag for "calibration runs" in which the number of photons per pixel is set to 0.
File size: 159.0 KB

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

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

Download in other formats: