Context Navigation

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

Visit:

Last change on this file since 2334 was 2334, checked in by blanch, 21 years ago
The concept of the camera prgoram has not changed but this version has quite a lot of changes to allow several Camera geometries as well as multitelescope. It is suposed to be the first working code for camera 0.7.
File size: 133.0 KB

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

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

Download in other formats: