Context Navigation

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

Visit:

Last change on this file since 5317 was 5305, checked in by moralejo, 20 years ago
Fixed error recently introduced in MLons::GetResponse. The NSB database was not correclt ycopied to the FADC when the randomly selected time was too close to the end of the database. This happened about 2% of times and produced some repetitive noise peaks in the FADC. Changed in camera.cxx the arguments of lons.SetSeed y lons_outer.SetSeed
File size: 152.6 KB

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

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

Download in other formats: