Context Navigation

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

Visit:

Last change on this file since 5407 was 5341, checked in by moralejo, 21 years ago
Removed unused variables fTelesTheta and fTelesPhi from MMcRunHeader. They were not useful because telescope orientation may change from event to event. Added fMirrorFraction to the MMcConfigRunHeader container in the camera output.
File size: 154.6 KB

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

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

Download in other formats: