Context Navigation

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

Visit:

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

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

Download in other formats: