Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

source: branches/start/MagicSoft/Simulation/Detector/TimeCam/timecam.cxx@ 13710

Visit:

Last change on this file since 13710 was 353, checked in by harald, 25 years ago
this is just the starting point of the further development of the camera project to simulate MonteCarloShower for MAGIC. The TimeCam is using also the information of the arrival time of the cerenkov photons (or the photoelectrons). This is the status of the program as presented on the general MAGIC meeting in Barcelona! Each one is invited to collaborate!!!
File size: 74.0 KB

Line
1	//!/////////////////////////////////////////////////////////////////////
2	//
3	// camera
4	//
5	// @file camera.cxx
6	// @title Camera simulation
7	// @subtitle Code for the simulation of the camera phase
8	// @desc Code for the simulation of the camera of CT1 and MAGIC
9	// @author J C Gonzalez
10	// @email gonzalez@mppmu.mpg.de
11	// @date Thu May 7 16:24:22 1998
12	//
13	//----------------------------------------------------------------------
14	//
15	// Created: Thu May 7 16:24:22 1998
16	// Author: Jose Carlos Gonzalez
17	// Purpose: Program for reflector simulation
18	// Notes: See files README for details
19	//
20	//----------------------------------------------------------------------
21	//
22	// $RCSfile: timecam.cxx,v $
23	// $Revision: 1.1.1.1 $
24	// $Author: harald $
25	// $Date: 2000-02-08 15:13:44 $
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 "TApplication.h"
52
53	#include "TFile.h"
54	#include "TTree.h"
55	#include "TBranch.h"
56	#include "TCanvas.h"
57
58	#include "MTrigger.hxx"
59	#include "MFadc.hxx"
60
61	#include "MRawEvt.h"
62	#include "MMcEvt.h"
63	#include "MMcTrig.hxx"
64
65	/*!@"
66
67	All the defines are located in the file \|camera.h\|.
68
69	@"*/
70
71	#include "timecam.h"
72	//!@}
73
74	/*!@"
75
76	The following set of flags are used in time of compilation. They do
77	not affect directly the behaviour of the program at run-time
78	(though, of course, if you disconnected the option for
79	implementation of the Trigger logic, you will not be able to use any
80	trigger at all. The 'default' values mean default in the sense of
81	what you got from the server when you obtained this program.
82
83	@"*/
84
85	//!@{
86
87	// flag for debugging (default: OFF )
88	#define __DEBUG__
89	#undef __DEBUG__
90
91	// flag for NNT in CT1 camera (default: ON )
92	#undef __CT1_NO_NEIGHBOURS__
93	#define __CT1_NO_NEIGHBOURS__
94
95	// flag for calculation of NSB (default: ON )
96	#undef __NSB__
97	#define __NSB__
98
99	// flag for calculation of QE for pixels (default: ON )
100	#undef __QE__
101	#define __QE__
102
103
104	// flag for implementation of DETAIL_TRIGGER (default: ON )
105	//
106	// This is the new implementation of Trigger studies
107	// It relies on a better simulation of the time stucture
108	// of the PhotoMultiplier. For more details see the
109	// documentation of the --> class MTrigger <--
110	#undef __DETAIL_TRIGGER__
111	#define __DETAIL_TRIGGER__
112
113	// flag for implementation of TRIGGER (default: ON )
114	#define __TRIGGER__
115	#undef __TRIGGER__
116
117	// flag for implementation of Tail-Cut (default: ON )
118	#define __TAILCUT__
119	#undef __TAILCUT__
120
121
122	// flag for calculation of islands stat. (default: ON )
123	#define __ISLANDS__
124	#undef __ISLANDS__
125
126	// flag for calculation of image parameters (default: ON )
127	#define __MOMENTS__
128	#undef __MOMENTS__
129
130	// flag for ROOT (default: ON )
131	#undef __ROOT__
132	#define __ROOT__
133
134	// flag for INTERAKTIV (default: OFF )
135	#undef __INTERAKTIV__
136	#define __INTERAKTIV__
137	//!@}
138
139	//=-----------------------------------------------------------
140	//!@subsection Definition of global variables.
141
142	/*!@"
143
144	Now we define some global variables with data about the telescope,
145	such as "focal distance", number of pixels/mirrors,
146	"size of the camera", and so on.
147
148	@"*/
149
150	/*!@"
151
152	Depending on the telescope we are using (CT1 or MAGIC), the
153	information stored in the definition file is different.
154	The variable \|ct_Type\| has the value 0 when we use
155	CT1, and 1 when we use MAGIC.
156
157	@"*/
158
159	//!@{
160	static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC
161	//!@}
162
163	/*!@"
164
165	And this is the information about the whole telescope.
166
167	@"*/
168
169	//!@{
170
171	// parameters of the CT (from the CT definition file)
172
173	////@: Focal distances [cm]
174	//static float *ct_Focal;
175
176	//@: Mean Focal distances [cm]
177	static float ct_Focal_mean;
178
179	//@: STDev. Focal distances [cm]
180	static float ct_Focal_std;
181
182	//@: Mean Point Spread function [cm]
183	static float ct_PSpread_mean;
184
185	//@: STDev. Point Spread function [cm]
186	static float ct_PSpread_std;
187
188	//@: STDev. Adjustmente deviation [cm]
189	static float ct_Adjustment_std;
190
191	//@: Radius of the Black Spot in mirror [cm]
192	static float ct_BlackSpot_rad;
193
194	//@: Radius of one mirror [cm]
195	static float ct_RMirror;
196
197	//@: Camera width [cm]
198	static float ct_CameraWidth;
199
200	//@: Pixel width [cm]
201	static float ct_PixelWidth;
202
203	//@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)
204	static float ct_PixelWidth_corner_2_corner;
205
206	//@: ct_PixelWidth_corner_2_corner / 2
207	static float ct_PixelWidth_corner_2_corner_half;
208
209	//@: Number of mirrors
210	static int ct_NMirrors = 0;
211
212	//@: Number of pixels
213	static int ct_NPixels;
214
215	//@: Number of pixels
216	static int ct_NCentralPixels;
217
218	//@: Number of pixels
219	static int ct_NGapPixels;
220
221	//@: ct_Apot = ct_PixelWidth / 2
222	static float ct_Apot;
223
224	//@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth
225	static float ct_2Apot;
226
227	//@: name of the CT definition file to use
228	static char ct_filename[256];
229
230	//@: list of showers to be skipped
231	static int *Skip;
232
233	//@: number of showers to be skipped
234	static int nSkip=0;
235
236	//@: flag: TRUE: data come from STDIN; FALSE: from file
237	static int Data_From_STDIN = FALSE;
238
239	//@: flag: TRUE: write all images to output; FALSE: only triggered showers
240	static int Write_All_Images = FALSE;
241
242	static int Write_McEvt = TRUE;
243	static int Write_McTrig = FALSE;
244	static int Write_RawEvt = FALSE;
245
246	//@: flag: TRUE: selection on the energy
247	static int Select_Energy = TRUE;
248
249	//@: Lower edge of the selected energy range (in GeV)
250	static float Select_Energy_le = 0.0;
251
252	//@: Upper edge of the selected energy range (in GeV)
253	static float Select_Energy_ue = 100000.0;
254
255	//!@}
256
257	/*!@"
258
259	The following double-pointer is a 2-dimensional table with information
260	about each pixel. The routine read_pixels will generate
261	the information for filling it using igen_pixel_coordinates().
262
263	@"*/
264
265	//!@{
266	// Pointer to a tables/Arrays with information about the pixels
267	// and data stored on them with information about the pixels
268
269	//@: table for IJ sys.
270	static float pixels[PIX_ARRAY_SIDE][PIX_ARRAY_SIDE][4];
271
272	//@: coordinates x,y for each pixel
273	static float **pixary;
274
275	//@: indexes of pixels neighbours of a given one
276	static int **pixneig;
277
278	//@: number of neighbours a pixel have
279	static int *npixneig;
280
281	//@: contents of the pixels (ph.e.)
282	static float *fnpix;
283
284	//@: contents of the pixels (ph.e.) after cleanning
285	static float *fnpixclean;
286
287	//!@}
288
289	/*!@"
290
291	The following double-pointer is a 2-dimensional table with the
292	Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
293	of the wavelength @$\lambda@$. The routine \|read_pixels()\| will read
294	also this information from the file \|qe.dat\|.
295
296	@"*/
297
298	//!@{
299	// Pointer to a table with QE, number of datapoints, and wavelengths
300
301	//@: table of QE
302	static float ***QE;
303
304	//@: number of datapoints for the QE curve
305	static int pointsQE;
306
307	//@: table of QE
308	static float *QElambda;
309	//!@}
310
311	/*!@"
312
313	The following double-pointer is a 2-dimensional table with information
314	about each mirror in the dish. The routine \|read_ct_file()\| will read
315	this information from the CT definition file.
316
317	@"*/
318
319	//!@{
320	// Pointer to a table with the following info.:
321
322	static float **ct_data;
323
324	/*
325	* TYPE=0 (CT1)
326	* i s rho theta x y z thetan phin xn yn zn
327	*
328	* i : number of the mirror
329	* s : arc length [cm]
330	* rho : polar rho of the position of the center of the mirror [cm]
331	* theta : polar angle of the position of the center of the mirror [cm]
332	* x : x coordinate of the center of the mirror [cm]
333	* y : y coordinate of the center of the mirror [cm]
334	* z : z coordinate of the center of the mirror [cm]
335	* thetan : polar theta angle of the direction where the mirror points to
336	* phin : polar phi angle of the direction where the mirror points to
337	* xn : xn coordinate of the normal vector in the center (normalized)
338	* yn : yn coordinate of the normal vector in the center (normalized)
339	* zn : zn coordinate of the normal vector in the center (normalized)
340	*
341	* TYPE=1 (MAGIC)
342	* i f sx sy x y z thetan phin
343	*
344	* i : number of the mirror
345	* f : focal distance of that mirror
346	* sx : curvilinear coordinate of mirror's center in X[cm]
347	* sy : curvilinear coordinate of mirror's center in X[cm]
348	* x : x coordinate of the center of the mirror [cm]
349	* y : y coordinate of the center of the mirror [cm]
350	* z : z coordinate of the center of the mirror [cm]
351	* thetan : polar theta angle of the direction where the mirror points to
352	* phin : polar phi angle of the direction where the mirror points to
353	* xn : xn coordinate of the normal vector in the center (normalized)
354	* yn : yn coordinate of the normal vector in the center (normalized)
355	* zn : zn coordinate of the normal vector in the center (normalized)
356	*/
357	//!@}
358
359	/*!@"
360
361	We define a table into where random numbers will be stored.
362	The routines used for random number generation are provided by
363	\|RANLIB\| (taken from NETLIB, \|www.netlib.org\|), and by
364	the routine \|double drand48(void)\| (prototype defined in
365	\|stdlib.h\|) through the macro \|RandomNumber\| defined in
366	\|camera.h\|.
367
368	@"*/
369
370	//!@{
371	// table of random numbers
372
373	// (unused)
374	//static double RandomNumbers[500];
375	//!@}
376
377	/*!@"
378
379	The following is a variable to count the number of Cphotons
380	in the different steps of the simulation.
381	The definition is as follows:
382	@[
383	\mbox{CountCphotons}[ \mbox{FILTER} ] \equiv
384	\mbox{\it Number of photons after the filter} \mbox{FILTER}
385	@]
386	The filters are defined and can be found in the file \|camera.h\|.
387
388	@"*/
389
390	//!@{
391	// vector to count photons at any given step of the simulation
392
393	//static int CountCphotons[10];
394	//!@}
395
396	/*!@"
397
398	The following are the set of parameters calculated for each image.
399	The routines for their calculations are in \|moments.cxx\|.
400
401	@"*/
402
403	//!@{
404	// parameters of the images
405
406	//static Moments_Info *moments_ptr;
407	//static LenWid_Info *lenwid_ptr;
408
409	//static float *maxs;
410	//static int *nmaxs;
411	//static float length, width, dist, xdist, azw, miss, alpha, *conc;
412	//static float phiasym, asymx, asymy;
413	//static float charge, smax, maxtrigthr_phe;
414
415	//!@}
416
417	extern char FileName[];
418
419
420	//=-----------------------------------------------------------
421	// @subsection Main program.
422
423	//!@{
424
425	//++++++++++++++++++++++++++++++++++++++++
426	// MAIN PROGRAM
427	//----------------------------------------
428
429	int main(int argc, char **argv)
430	{
431	//!@' @#### Definition of variables.
432	//@'
433
434	char inname[256]; //@< input file name
435	char rootname[256] ; //@< ROOT file name
436
437	char parname[256]; //@< parameters file name
438
439	char sign[20]; //@< initialize sign
440
441	char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
442
443	ifstream inputfile; //@< stream for the input file
444
445	MCEventHeader mcevth; //@< Event Header class (MC)
446	MCCphoton cphoton; //@< Cherenkov Photon class (MC)
447
448	float thetaCT, phiCT; //@< parameters of a given shower
449	float thetashw, phishw; //@< parameters of a given shower
450	float coreD, coreX, coreY; //@< core position and distance
451	float impactD; //@< impact parameter
452	float l1, m1, n1; //@< auxiliary variables
453	float l2, m2, n2; //@< auxiliary variables
454	float num, den; //@< auxiliary variables
455
456	int nshow=0; //@< partial number of shower in a given run
457	int ntshow=0; //@< total number of showers
458	int ncph_in=0; //@< number of (input) photons in one shower
459	int ntcph_in=0; //@< total number of input photons
460
461	int ncph_out=0; //@< number of (output) photons in one shower
462	int ntcph_out=0; //@< total number of output photons
463
464	int i, ii, k; //@< simple counters
465
466	float t_ini; //@< time of the first Cphoton read in
467	float t; //@< time for a single photon
468	float t_first ; //@< time of the first cerenkov particle reaching the ground
469	float t_last ; //@< time of the last cerenkov particle reaching the ground
470
471	int t_chan ; //@< the bin (channel) in time of a single photon
472
473	int startchan ; //@< the first bin with entries in the time slices
474
475	float cx, cy, ci, cj; //@< coordinates in the XY and IJ systems
476	int ici, icj, iici, iicj; //@< coordinates in the IJ (integers)
477
478	int nPMT; //@< number of pixel
479
480	float wl, last_wl; //@< wavelength of the photon
481	float qe; //@< quantum efficiency
482	float **qeptr;
483
484	int simulateNSB; //@< Will we simulate NSB?
485	float meanNSB; //@< NSB mean value (per pixel and ns)
486	float meanPois; //@< NSB mean value for simulation
487	float timeNSB ; //@< time of NSB photon
488
489	int nIslandsCut; //@< Islands Cut value
490	int countIslands; //@< Will we count the islands?
491	int anaPixels;
492
493	float fCorrection; //@< Factor to apply to pixel values (def. 1.)
494
495	// int ntrigger = 0; //@< number of triggers in the whole file
496
497	float plateScale_cm2deg; //@< plate scale (deg/cm)
498	float degTriggerZone; //@< trigger area in the camera (radius, in deg.)
499
500	float dtheta, dphi; //@< deviations of CT from shower axis
501
502	int still_in_loop = FALSE;
503
504	char Signature[20];
505
506	float *image_data;
507	int nvar;
508	// int hidt;
509
510	struct camera cam; // structure holding the camera definition
511
512	//!@' @#### Definition of variables for \|getopt()\|.
513	//@'
514
515	int ch, errflg = 0; //@< used by getopt
516
517	/*!@'
518
519	@#### Beginning of the program.
520
521	We start with the main program. First we (could) make some
522	presentation, and follows the reading of the parameters file (now
523	from the \|stdin\|), the reading of the CT parameters file, and the
524	creation of the output file, where the processed data will be
525	stored.
526
527	*/
528
529	//++
530	// START
531	//--
532
533	// make unbuffered output
534
535	cout.setf ( ios::stdio );
536
537	// parse command line options (see reflector.h)
538
539	parname[0] = '\0';
540
541	optarg = NULL;
542	while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
543	switch (ch) {
544	case 'f':
545	strcpy(parname, optarg);
546	break;
547	case 'h':
548	usage();
549	break;
550	default :
551	errflg++;
552	}
553
554	// show help if error
555
556	if ( errflg>0 )
557	usage();
558
559	// make some sort of presentation
560
561	present();
562
563	// read parameters file
564
565	if ( strlen(parname) < 1 )
566	readparam(NULL);
567	else
568	readparam(parname);
569
570	// read data from file or from STDIN?
571
572	Data_From_STDIN = get_data_from_stdin();
573
574	// write all images, even those without trigger?
575
576	Write_All_Images = get_write_all_images();
577
578	Write_McEvt = get_write_McEvt() ;
579	Write_McTrig = get_write_McTrig() ;
580	Write_RawEvt = get_write_RawEvt() ;
581
582	// get filenames
583
584	strcpy( inname, get_input_filename() );
585	strcpy( rootname, get_root_filename() );
586	strcpy( ct_filename, get_ct_filename() );
587
588	// get different parameters of the simulation
589
590	simulateNSB = get_nsb( &meanNSB );
591	countIslands = get_islands_cut( &nIslandsCut );
592
593	// get selections on the parameters
594
595	Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
596
597	// log filenames information
598
599	log(SIGNATURE,
600	"%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n",
601	"Filenames",
602	"In", inname,
603	"ROOT", rootname,
604	"CT", ct_filename);
605
606
607	// log flags information
608
609	log(SIGNATURE,
610	"%s:\n\t%20s: %s\n\t%20s: %s\n",
611	"Flags",
612	"Data_From_STDIN", ONoff(Data_From_STDIN),
613	"Write_All_Images", ONoff(Write_All_Images));
614
615
616	// log flags information
617
618	log(SIGNATURE,
619	"%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
620	"Root output",
621	"Write_McEvt", ONoff(Write_McEvt),
622	"Write_McTrig", ONoff(Write_McTrig),
623	"Write_RawEvt", ONoff(Write_RawEvt));
624
625	// log parameters information
626
627	log(SIGNATURE,
628	"%s:\n\t%20s: %f %s\n\t%20s: %f %s\n",
629	"Parameters",
630	"NSB (phes/pixel)", meanNSB, ONoff(simulateNSB),
631	"i0 (Islands-cut)", nIslandsCut, ONoff(countIslands));
632
633	// log selections
634
635	log(SIGNATURE,
636	"%s:\n\t%20s: %s (%f:%f)\n",
637	"Selections:",
638	"Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
639
640	// set all random numbers seeds
641
642	setall( get_seeds(0), get_seeds(1) );
643
644	// get list of showers to evt. skip
645
646	nSkip = get_nskip_showers();
647
648	if (nSkip > 0) {
649	Skip = new int[ nSkip ];
650	get_skip_showers( Skip );
651
652	log(SIGNATURE, "There are some showers to skip:\n");
653	for (i=0; i<nSkip; ++i)
654	log(SIGNATURE, "\tshower # %d\n", Skip[i]);
655	}
656
657	// read parameters from the ct.def file
658
659	read_ct_file();
660
661	// read pixels data
662
663	read_pixels(&cam);
664
665
666	Int_t Lev0, Lev1, Lev2 ;
667
668	// initialise ROOT
669
670	TROOT simple("simple", "MAGIC Telescope Monte Carlo");
671
672
673	#ifdef __DETAIL_TRIGGER__
674
675	MTrigger Trigger ; //@< A instance of the Class MTrigger
676
677	MMcTrig *McTrig = new MMcTrig() ;
678
679
680	// MFadc fadc ;
681
682	#endif // __DETAIL_TRIGGER__
683
684	#ifdef __ROOT__
685
686	MRawEvt *Evt = new MRawEvt() ;
687	MMcEvt *McEvt = new MMcEvt ();
688
689	// initalize the ROOT file
690	//
691	// erzeuge ein Root file
692	//
693
694	TFile outfile ( rootname , "RECREATE" ) ;
695
696	//
697	// create a Tree for the Event data stream
698	//
699
700	TTree EvtTree("EvtTree","Events of Run");
701
702	Int_t bsize=128000; Int_t split=1;
703
704	//
705	// check which branches to create (you are able to controll this via
706	// the input file
707	//
708
709	if ( Write_McEvt == TRUE ) {
710	EvtTree.Branch("MMcEvt","MMcEvt",
711	&McEvt, bsize, split);
712	}
713
714	if ( Write_McTrig == TRUE ) {
715	EvtTree.Branch("MMcTrig","MMcTrig",
716	&McTrig, bsize, split);
717	}
718
719	if ( Write_RawEvt == TRUE ) {
720	EvtTree.Branch("MRawEvt","MRawEvt",
721	&Evt, bsize, split);
722	}
723
724
725	unsigned short ulli = 0 ;
726
727	#endif // __ROOT__
728
729	#ifdef __NSB__
730	TRandom GenNSB ;
731	#endif // __NSB__
732
733	#ifdef __INTERAKTIV__
734	TApplication theApp("App", &argc, argv);
735
736	if (gROOT->IsBatch()) {
737	fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
738	// return 1;
739	}
740	#endif // __INTERAKTIV__
741
742
743	// for safety and for dimensioning image_data: count the elements in the
744	// diagnostic data branch
745
746	i=0;
747	i++; // "n"
748	i++; // "primary"
749	i++; // "energy"
750	i++; // "cored"
751	i++; // "impact"
752	i++; // "xcore"
753	i++; // "ycore"
754	i++; // "theta"
755	i++; // "phi"
756	i++; // "deviations"
757	i++; // "dtheta"
758	i++; // "dphi"
759	i++; // "trigger"
760	i++; // "ncphs"
761	i++; // "maxpassthr_phe"
762	i++; // "nphes"
763	i++; // "nphes2"
764	i++; // "length"
765	i++; // "width"
766	i++; // "dist"
767	i++; // "xdist"
768	i++; // "azw"
769	i++; // "miss"
770	i++; // "alpha"
771	i++; // "conc2"
772	i++; // "conc3"
773	i++; // "conc4"
774	i++; // "conc5"
775	i++; // "conc6"
776	i++; // "conc7"
777	i++; // "conc8"
778	i++; // "conc9"
779	i++; // "conc10"
780	i++; // "asymx"
781	i++; // "asymy"
782	i++; // "phiasym"
783
784	nvar = i;
785	image_data = new float[nvar];
786
787	// set plate scale (deg/cm) and trigger area (deg)
788
789	plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;
790
791	if ( ! get_trigger_radius( &degTriggerZone ) )
792	degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
793
794	if ( ! get_correction( &fCorrection ) )
795	fCorrection = 1.0;
796
797	// number of pixels for parameters
798
799	anaPixels = get_ana_pixels();
800	anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
801
802	// open input file if we DO read data from a file
803
804	if (! Data_From_STDIN) {
805	log( SIGNATURE, "Openning input \"rfl\" file %s\n", inname );
806	inputfile.open( inname );
807	if ( inputfile.bad() )
808	error( SIGNATURE, "Cannot open input file: %s\n", inname );
809	}
810
811	// get signature, and check it
812	// NOTE: this part repeats further down in the code;
813	// if you change something here you probably want to change it
814	// there as well
815
816	strcpy(Signature, REFL_SIGNATURE);
817
818	strcpy(sign, Signature);
819
820	if ( Data_From_STDIN )
821	cin.read( (char *)sign, strlen(Signature));
822	else
823	inputfile.read( (char *)sign, strlen(Signature));
824
825	if (strcmp(sign, Signature) != 0) {
826	cerr << "ERROR: Signature of .rfl file is not correct\n";
827	cerr << '"' << sign << '"' << '\n';
828	cerr << "should be: " << Signature << '\n';
829	exit(1);
830	}
831
832	if ( Data_From_STDIN )
833	cin.read( (char *)sign, 1);
834	else
835	inputfile.read( (char *)sign, 1);
836
837	// initializes flag
838
839	strcpy( flag, " \0" );
840
841	// allocate space for PMTs numbers of pixels
842
843	fnpix = new float [ ct_NPixels ];
844	fnpixclean = new float [ ct_NPixels ];
845
846	#ifdef __ROOT__
847
848	// construct MFADC
849
850	#endif // __ROOT__
851
852
853	// moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
854
855	//!@' @#### Main loop.
856	//@'
857
858	//begin my version
859
860	// get flag
861
862	if ( Data_From_STDIN )
863	cin.read( flag, SIZE_OF_FLAGS );
864	else
865	inputfile.read ( flag, SIZE_OF_FLAGS );
866
867	// loop over the file
868
869	still_in_loop = TRUE;
870
871	while (
872	((! Data_From_STDIN) && (! inputfile.eof()))
873	\|\|
874	(Data_From_STDIN && still_in_loop)
875	) {
876
877
878
879	// reading .rfl files
880	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
881	log(SIGNATURE, "End of file . . .\n");
882	still_in_loop = FALSE;
883
884	if ( Data_From_STDIN )
885	cin.read( (char *)sign, 1);
886	else
887	inputfile.read( (char *)sign, 1);
888
889
890	}
891	else if(!isA( flag, FLAG_START_OF_RUN )){
892	error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
893	}
894	else { // found start of run
895	nshow=0;
896	// read next flag
897
898	if ( Data_From_STDIN )
899	cin.read( flag, SIZE_OF_FLAGS );
900	else
901	inputfile.read ( flag, SIZE_OF_FLAGS );
902
903	while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
904	/*!@'
905
906	For the case \|FLAG\_START\_OF\_EVENT\|,
907	we read each Cherenkov photon, and follow these steps:
908
909	@enumerate
910
911	@- Transform XY-coordinates to IJ-coordinates.
912
913	@- With this, we obtain the pixel where the photon hits.
914
915	@- Use the wavelength $\lambda$ and the table of QE, and
916	calculate the estimated (third order interpolated) quantum
917	efficiency for that photon. The photon can be rejected.
918
919	@- If accepted, then add to the pixel.
920
921	@endenumerate
922
923	In principle, we should stop here, and use another program to
924	'smear' the image, to add the Night Sky Background, and to
925	simulate the trigger logic, but we will make this program
926	quick and dirty, and include all here.
927
928	If we are reading PHE files, we jump to the point where the
929	pixelization process already has finished.
930
931	*/
932
933	++nshow;
934
935	if ( fmod ( nshow, 1000. ) == 0. )
936	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
937
938
939	// get MCEventHeader
940
941	if ( Data_From_STDIN )
942	cin.read( (char*)&mcevth, mcevth.mysize() );
943	else
944	mcevth.read( inputfile );
945
946	// calculate core distance and impact parameter
947
948	coreD = mcevth.get_core(&coreX, &coreY);
949
950	// calculate impact parameter (shortest distance betwee the original
951	// trajectory of the primary (assumed shower-axis) and the
952	// direction where the telescope points to
953	//
954	// we use the following equation, given that the shower core position
955	// is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
956	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
957	// and (l2,m2,n2)
958	//
959	// \| \|
960	// \| x1-x2 y1-y2 z1-z2 \|
961	// \| \|
962	// + \| l1 m1 n1 \|
963	// - \| \|
964	// \| l2 m2 n2 \|
965	// \| \|
966	// dist = ------------------------------------ ( > 0 )
967	// [ \|l1 m1\|2 \|m1 n1\|2 \|n1 l1\|2 ]1/2
968	// [ \| \| + \| \| + \| \| ]
969	// [ \|l2 m2\| \|m2 n2\| \|n2 l2\| ]
970	//
971	// playing a little bit, we get this reduced for in our case:
972	//
973	//
974	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
975	// [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
976	// 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
977
978	// read the direction of the incoming shower
979
980	thetashw = mcevth.get_theta();
981	phishw = mcevth.get_phi();
982
983	// calculate vector for shower
984
985	l1 = sin(thetashw)*cos(phishw);
986	m1 = sin(thetashw)*sin(phishw);
987	n1 = cos(thetashw);
988
989	// read the deviation of the telescope with respect to the shower
990
991	mcevth.get_deviations ( &thetaCT, &phiCT );
992
993	if ( (thetaCT == 0.) && (phiCT == 0.) ) {
994
995	// CT was looking to the source (both lines are parallel)
996	// therefore, we calculate the impact parameter as the distance
997	// between the CT axis and the core position
998
999	impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1000
1001	} else {
1002
1003	// the shower comes off-axis
1004
1005	// obtain with this the final direction of the CT
1006
1007	thetaCT += thetashw;
1008	phiCT += phishw;
1009
1010	// calculate vector for telescope
1011
1012	l2 = sin(thetaCT)*cos(phiCT);
1013	m2 = sin(thetaCT)*sin(phiCT);
1014	n2 = cos(thetaCT);
1015
1016	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1017	den = (SQR(l1m2 - l2m1) +
1018	SQR(m1n2 - m2n1) +
1019	SQR(n1l2 - n2l1));
1020	den = sqrt(den);
1021
1022	impactD = fabs(num)/den;
1023
1024	// fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ",
1025	// thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2,
1026	// num, den);
1027
1028	}
1029
1030	// clear camera
1031
1032	for ( i=0; i<ct_NPixels; ++i ){
1033
1034	fnpix[i] = 0.0;
1035	}
1036
1037	#ifdef __ROOT__
1038	// clear the MFADC
1039
1040	#endif // __ROOT__
1041
1042	ntcph_in +=ncph_in;
1043	ncph_in = 0;
1044
1045	ntcph_out +=ncph_out;
1046	ncph_out = 0;
1047
1048
1049	#ifdef __DETAIL_TRIGGER__
1050	//
1051	// clear Trigger
1052	//
1053
1054	Trigger.Reset() ;
1055	#endif // __DETAIL_TRIGGER__
1056
1057	// fadc.Reset() ;
1058
1059	//
1060	// Read out the first and last time of cerenkovs
1061	//
1062
1063	mcevth.get_times ( &t_first, &t_last ) ;
1064
1065	//- - - - - - - - - - - - - - - - - - - - - - - - -
1066	// read photons and "map" them into the pixels
1067	//--------------------------------------------------
1068
1069	// initialize CPhoton
1070
1071	cphoton.fill(0., 0., 0., 0., 0., 0., 0., 0.);
1072
1073	// read the photons data
1074
1075	if ( Data_From_STDIN )
1076	cin.read( flag, SIZE_OF_FLAGS );
1077	else
1078	inputfile.read ( flag, SIZE_OF_FLAGS );
1079
1080	// loop over the photons
1081
1082	t_ini = -99999;
1083
1084	while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
1085
1086	memcpy( (char*)&cphoton, flag, SIZE_OF_FLAGS );
1087
1088	if ( Data_From_STDIN )
1089	cin.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
1090	else
1091	inputfile.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
1092
1093	// increase number of photons
1094
1095	ncph_in++;
1096
1097	t = cphoton.get_t() ;
1098
1099
1100	/*!@'
1101
1102	@#### Pixelization (for the central pixels).
1103
1104	In order to calculate the coordinates, we use the
1105	change of system described in the documentation
1106	of the source code of \|pixel\_coord.cxx\|.
1107	Then, we will use simply the matrix of change
1108	from one system to the other. In our case, this is:
1109
1110	@[
1111	\begin{bmatrix}X\\Y\\\end{bmatrix}
1112	=
1113	\begin{bmatrix}
1114	1 & \cos(60^\circ)\\
1115	0 & \sin(60^\circ)\\
1116	\end{bmatrix}
1117	\begin{bmatrix}I\\J\\\end{bmatrix}
1118	@]
1119
1120	and hence
1121
1122	@[
1123	\begin{bmatrix}I\\J\\\end{bmatrix}
1124	=
1125	\begin{bmatrix}
1126	1 & -\frac{\cos(60^\circ)}{\sin(60^\circ)}\\
1127	0 &\frac{1}{\sin(60^\circ)}\\
1128	\end{bmatrix}
1129	\begin{bmatrix}X\\Y\\\end{bmatrix}
1130	@]
1131
1132	*/
1133
1134	//+++
1135	// Pixelization
1136	//---
1137
1138	// calculate ij-coordinates
1139
1140	// We use a change of coordinate system, using the following
1141	// matrix of change (m^-1) (this is taken from Mathematica output).
1142	/*
1143	* In[1]:= m={{1,cos60},{0,sin60}}; MatrixForm[m]
1144	*
1145	* Out[1]//MatrixForm= 1 cos60
1146	*
1147	* 0 sin60
1148	*
1149	* In[2]:= inv=Inverse[m]; MatrixForm[inv]
1150	*
1151	* Out[2]//MatrixForm= cos60
1152	* -(-----)
1153	* 1 sin60
1154	*
1155	* 1
1156	* -----
1157	* 0 sin60
1158	*
1159	*/
1160
1161	// go to IJ-coordinate system
1162
1163	cx = cphoton.get_x();
1164	cy = cphoton.get_y();
1165
1166	// get wavelength
1167
1168	last_wl = wl;
1169	wl = cphoton.get_wl();
1170
1171	if ( wl < 1.0 )
1172	break;
1173
1174	if ( (wl > 600.0) \|\| (wl < 290.0) )
1175	break;
1176
1177	// check if photon is inside outermost camera radius
1178
1179	if(sqrt(cxcx + cycy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)){
1180
1181	// read next CPhoton
1182	if ( Data_From_STDIN )
1183	cin.read( flag, SIZE_OF_FLAGS );
1184	else
1185	inputfile.read ( flag, SIZE_OF_FLAGS );
1186
1187	// go to beginning of loop, the photon is lost
1188	continue;
1189
1190	}
1191
1192	// cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
1193
1194	ci = floor( (cx - cy*COS60/SIN60)/ ct_2Apot + 0.5);
1195	cj = floor( (cy/SIN60) / ct_2Apot + 0.5);
1196
1197	ici = (int)(ci);
1198	icj = (int)(cj);
1199
1200	iici = ici+PIX_ARRAY_HALF_SIDE;
1201	iicj = icj+PIX_ARRAY_HALF_SIDE;
1202
1203	// is it inside the array?
1204
1205	if ( (iici > 0) && (iici < PIX_ARRAY_SIDE) &&
1206	(iicj > 0) && (iicj < PIX_ARRAY_SIDE) ) {
1207
1208	// try to put into pixel
1209
1210	// obtain the pixel number for this photon
1211
1212	nPMT = (int)
1213	pixels[ici+PIX_ARRAY_HALF_SIDE][icj+PIX_ARRAY_HALF_SIDE][PIXNUM];
1214
1215	}
1216	else{
1217
1218	nPMT = -1;
1219
1220	}
1221
1222	// check if outside the central camera
1223
1224	if ( (nPMT < 0) \|\| (nPMT >= ct_NCentralPixels) ) {
1225
1226	// check the outer pixels
1227	nPMT = -1;
1228
1229	for(i=ct_NCentralPixels; i<ct_NPixels; i++){
1230	if( bpoint_is_in_pix( cx, cy, i, &cam) ){
1231	nPMT = i;
1232	break;
1233	}
1234	}
1235
1236	if(nPMT==-1){// the photon is in none of the pixels
1237
1238	// read next CPhoton
1239	if ( Data_From_STDIN )
1240	cin.read( flag, SIZE_OF_FLAGS );
1241	else
1242	inputfile.read ( flag, SIZE_OF_FLAGS );
1243
1244	// go to beginning of loop, the photon is lost
1245	continue;
1246	}
1247
1248	}
1249
1250	#ifdef __QE__
1251
1252	//!@' @#### QE simulation.
1253	//@'
1254
1255	//+++
1256	// QE simulation
1257	//---
1258
1259	// find data point to be used in Lagrange interpolation (-> k)
1260
1261	qeptr = (float **)QE[nPMT];
1262
1263	FindLagrange(qeptr,k,wl);
1264
1265	// if random > quantum efficiency, reject it
1266
1267	qe = Lagrange(qeptr,k,wl) / 100.0;
1268
1269	// fprintf(stdout, "%f\n", qe);
1270
1271	if ( RandomNumber > qe ) {
1272
1273	// read next CPhoton
1274	if ( Data_From_STDIN )
1275	cin.read( flag, SIZE_OF_FLAGS );
1276	else
1277	inputfile.read ( flag, SIZE_OF_FLAGS );
1278
1279	// go to beginning of loop
1280	continue;
1281
1282	}
1283
1284	#endif // __QE__
1285
1286	//+++
1287	// Cphoton is accepted
1288	//---
1289
1290	ncph_out++ ;
1291
1292	// increase the number of Cphs. in the PMT, i.e.,
1293	// increase in one unit the counter of the photons
1294	// stored in the pixel nPMT
1295
1296	fnpix[nPMT] += 1.0;
1297
1298	#ifdef __DETAIL_TRIGGER__
1299	//
1300	// fill the Trigger class with this phe
1301	//
1302	//
1303	Trigger.Fill( nPMT, ( t - t_first ) ) ;
1304
1305	// fadc.Fill( nPMT, ( t - t_first ), Trigger.Fill( nPMT, ( t - t_first ) ) ) ;
1306
1307	#endif // __DETAIL_TRIGGER__
1308
1309	// read next CPhoton
1310
1311	if ( Data_From_STDIN )
1312	cin.read( flag, SIZE_OF_FLAGS );
1313	else
1314	inputfile.read ( flag, SIZE_OF_FLAGS );
1315
1316	} // end while, i.e. found end of event
1317
1318	if ( fmod ( nshow, 1000. ) == 0. )
1319	log(SIGNATURE,
1320	"End of this event: in: %d cphs(+%d). out: %d cphs(+%d). .\n",
1321	ncph_in, ntcph_in,
1322	ncph_out, ntcph_out);
1323
1324	// show number of photons
1325
1326	//cout << ncph_in << " photons read . . . " << endl << flush;
1327
1328	// skip it ?
1329
1330	for ( i=0; i<nSkip; ++i ) {
1331	if (Skip[i] == (nshow+ntshow)) {
1332	i = -1;
1333	break;
1334	}
1335	}
1336
1337	// if after the previous loop, the exit value of i is -1
1338	// then the shower number is in the list of showers to be
1339	// skipped
1340
1341	if (i == -1) {
1342	log(SIGNATURE, "\t\tskipped!\n");
1343	continue;
1344	}
1345
1346	/*!@'
1347
1348	After reading all the Cherenkov photons for a given event,
1349	we have in the table of number of photons for each pixel
1350	only the 'raw' amount of Cherenkov photons @$n_p@$. Now, we
1351	should take this number as the mean value of the
1352	distribution of photons in that pixel @$p@$, following a
1353	Poisson distribution.
1354
1355	@[ n_p \equiv \mu_p @]
1356
1357	and with this number the amount of light coming from the
1358	shower is calculated @$\hat{n}_p@$.
1359
1360	Then, we calculate the amount of Night Sky Background we
1361	must introduce in that pixel @$p@$. We calculate this using
1362	again a Poisson distribution with mean @$\mu_\mathrm{NSB}@$
1363	(defined in the \|camera.h\| file). The value of
1364	@$\mu_\mathrm{NSB}@$ is obtained from measurements. With this
1365	value, the amount of photons @$\hat{n}_\mathrm{NSB}@$ coming
1366	from the Night Sky Background is calculated.
1367
1368	Finally, the amount of photons for that pixels is:
1369	@[ \hat{n}_p^\mathrm{final} = \hat{n}_p + \hat{n}_\mathrm{NSB} @]
1370
1371	*/
1372
1373	// after reading all the photons, our camera is filled
1374
1375	if ( Select_Energy ) {
1376	if (( mcevth.get_energy() < Select_Energy_le ) \|\|
1377	( mcevth.get_energy() > Select_Energy_ue )) {
1378	log(SIGNATURE, "select_energy: shower rejected.\n");
1379	continue;
1380	}
1381	}
1382
1383	#ifdef __NSB__
1384
1385	//!@' @#### NSB (Night Sky Background) simulation.
1386	//@'
1387
1388	//+++
1389	// NSB simulation
1390	//---
1391
1392	// add NSB "noise"
1393	// TO DO: make meanNSB an array and read the contents from a file!
1394
1395	// if ( simulateNSB )
1396	// for ( i=0; i<ct_NPixels; ++i )
1397	// fnpix[i] += (float)ignpoi( meanNSB );
1398	// old version of Jose Carlos
1399
1400	if ( simulateNSB) {
1401	//
1402	// loop over all pixels and scramble the number
1403	// of NSB photons to put in it. For the number use
1404	// a poison distribution with a mean calculated from meanNSB
1405	// and the TOTAL_TRIGGER_TIME
1406	//
1407	for ( Int_t nsbPix = 0 ; nsbPix < CAMERA_PIXELS ; nsbPix++ ) {
1408	//
1409	//
1410	meanPois = meanNSB * TOTAL_TRIGGER_TIME ;
1411
1412	// loop over the scrambled number if Photons in this pixels
1413
1414	for ( Int_t photNSB=0; photNSB<GenNSB.Poisson(meanPois);
1415	photNSB++){
1416	//
1417	// now scramble the time at that the photo electron of the
1418	// NSB photon is leaving the photo cathod
1419	//
1420
1421	timeNSB = GenNSB.Rndm() * TOTAL_TRIGGER_TIME ;
1422
1423	Trigger.FillNSB ( nsbPix, timeNSB ) ;
1424
1425	}
1426
1427
1428	}
1429
1430	}
1431	#endif // __NSB__
1432
1433	// if we should apply any kind of correction, do it here.
1434
1435	for ( i=0; i<ct_NPixels; ++i )
1436	fnpix[i] *= fCorrection;
1437
1438	#ifdef __DETAIL_TRIGGER__
1439	//
1440	// now the noise of the electronic
1441	// (preamps, optical transmission,..) is introduced.
1442	// This is done inside the class MTrigger by the method ElecNoise.
1443	//
1444	Trigger.ElecNoise() ;
1445
1446	//
1447	// look if in all the signals in the trigger signal branch
1448	// is a possible Trigger. Therefore we habe to diskriminate all
1449	// the simulated analog signals (Method Diskriminate in class
1450	// MTrigger). We look simultanously for the moments at which
1451	// there are more than TRIGGER_MULTI pixels above the
1452	// CHANNEL_THRESHOLD.
1453	//
1454
1455	McTrig->SetZeroLevel( Lev0 = (Short_t) Trigger.Diskriminate() ) ;
1456	Lev1 = Lev2 = 0 ;
1457
1458	//
1459	// Start the First Level Trigger simulation
1460	//
1461
1462	if ( Lev0 > 0 ) {
1463	McTrig->SetFirstLevel ( Lev1 = Trigger.FirstLevel() ) ;
1464	}
1465
1466	#endif // __DETAIL_TRIGGER__
1467
1468	#ifdef __ROOT__
1469
1470	//
1471	// Fill the header of this event
1472	//
1473
1474	Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ;
1475
1476	//
1477	// fill the MMcEvt with all information
1478	//
1479
1480	McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1481	mcevth.get_energy(),
1482	mcevth.get_theta(),
1483	mcevth.get_phi(),
1484	mcevth.get_core(),
1485	mcevth.get_coreX(),
1486	mcevth.get_coreY(),
1487	impactD,
1488	ulli, ulli,
1489	(UShort_t) ncph_in,
1490	ulli,
1491	(UShort_t) ncph_out ) ;
1492
1493	//
1494	// write it out to the file outfile
1495	//
1496
1497	EvtTree.Fill() ;
1498
1499	#endif // __ROOT__
1500
1501	//
1502	// if a first level trigger occurred, then
1503	// 1. do some other stuff (not implemented)
1504	// 2. start the gui tool
1505
1506	if ( Lev1 > 0 ) {
1507
1508	// fadc.Scan( Trigger.GetFirstLevelTime(0) ) ;
1509
1510	#ifdef __INTERAKTIV__
1511	Trigger.ShowSignal(McEvt) ;
1512	#endif
1513	}
1514
1515
1516
1517	#ifdef __ROOT__
1518	// clear all
1519	Evt->Clear() ;
1520	McEvt->Clear() ;
1521	McTrig->Clear() ;
1522	#endif // __ROOT__
1523
1524
1525	//++++++++++++++++++++++++++++++++++++++++++++++++++
1526	// at this point we have a camera full of
1527	// ph.e.s
1528	// we should first apply the trigger condition,
1529	// and if there's trigger, then clean the image,
1530	// calculate the islands statistics and the
1531	// other parameters of the image (Hillas' parameters
1532	// and so on).
1533	//--------------------------------------------------
1534
1535	#ifdef __DEBUG__
1536	printf("\n");
1537
1538	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1539
1540	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1541
1542	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1543
1544	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1545
1546	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1547	(int)pixels[ici][icj][PIXNUM],
1548	pixels[ici][icj][PIXX],
1549	pixels[ici][icj][PIXY],
1550	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1551
1552	}
1553	}
1554	}
1555	}
1556
1557	for (i=0; i<ct_NPixels; ++i) {
1558	printf("%d (%d): ", i, npixneig[i]);
1559	for (j=0; j<npixneig[i]; ++i)
1560	printf(" %d", pixneig[i][j]);
1561	printf("\n");
1562	}
1563
1564	#endif // __DEBUG__
1565
1566
1567	//!@' @#### Save data.
1568	//@'
1569
1570	//++++++++++++++++++++++++++++++++++++++++++++++++++
1571	// we now have all information we want
1572	// the only thing we must do now is writing it to
1573	// the output file
1574	//--------------------------------------------------
1575
1576	//++
1577	// save the image to the file
1578	//--
1579
1580	if ( Data_From_STDIN )
1581	cin.read( flag, SIZE_OF_FLAGS );
1582	else
1583	inputfile.read ( flag, SIZE_OF_FLAGS );
1584
1585	} // end while there is a next event
1586
1587	if( !isA( flag, FLAG_END_OF_RUN )){
1588	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1589	}
1590	else { // found end of run
1591	ntshow += nshow;
1592	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1593
1594	// if ( Data_From_STDIN )
1595	// cin.read( flag, SIZE_OF_FLAGS );
1596	// else
1597	// inputfile.read ( flag, SIZE_OF_FLAGS );
1598
1599	// huschel start here
1600
1601	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1602	log(SIGNATURE, "End of file . . .\n");
1603	still_in_loop = FALSE;
1604
1605	if ((! Data_From_STDIN) && (! inputfile.eof())){
1606
1607	// we have concatenated input files.
1608	// get signature of the next part and check it.
1609	// NOTE: this part repeats further up in the code;
1610	// if you change something here you probably want to change it
1611	// there as well
1612
1613	strcpy(Signature, REFL_SIGNATURE);
1614
1615	strcpy(sign, Signature);
1616
1617	inputfile.read( (char *)sign, strlen(Signature));
1618
1619	if (strcmp(sign, Signature) != 0) {
1620	cerr << "ERROR: Signature of .rfl file is not correct\n";
1621	cerr << '"' << sign << '"' << '\n';
1622	cerr << "should be: " << Signature << '\n';
1623	exit(1);
1624	}
1625
1626	if ( Data_From_STDIN )
1627	cin.read( (char *)sign, 1);
1628	else
1629	inputfile.read( (char *)sign, 1);
1630
1631	}
1632
1633	// huschel ends here
1634
1635	} // end if found end of file
1636
1637	} // end if found end of run
1638
1639	if ( Data_From_STDIN )
1640	cin.read( flag, SIZE_OF_FLAGS );
1641	else
1642	inputfile.read ( flag, SIZE_OF_FLAGS );
1643
1644	} // end if else found start of run
1645	} // end big while loop
1646
1647	//!@' @#### End of program.
1648	//@'
1649
1650	//end my version
1651
1652	#ifdef __ROOT__
1653	//++
1654	// put the Event to the root file
1655	//--
1656
1657	EvtTree.Write() ;
1658	outfile.Write() ;
1659	outfile.Close() ;
1660
1661	#endif // __ROOT__
1662
1663	// close input file
1664
1665	ntcph_in += ncph_in;
1666	ntcph_out += ncph_out;
1667	log( SIGNATURE,
1668	"%d event(s), with a total of %d C.photons in and %d C.photons out \n",
1669	ntshow, ntcph_in, ntcph_out );
1670
1671	// log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
1672	// ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
1673
1674	// close files
1675
1676	log( SIGNATURE, "Closing files\n" );
1677
1678	inputfile.close();
1679
1680	#ifdef __DETAIL_TRIGGER__
1681	// Output of Trigger statistics
1682	//
1683
1684	// Trigger.PrintStat() ;
1685	#endif // __DETAIL_TRIGGER__
1686
1687	// program finished
1688
1689	log( SIGNATURE, "Done.\n");
1690
1691	return( 0 );
1692
1693	}
1694	//!@}
1695
1696	// @T \newpage
1697
1698	//!@subsection Functions definition.
1699
1700	//!-----------------------------------------------------------
1701	// @name present
1702	//
1703	// @desc Make some presentation
1704	//
1705	// @date Sat Jun 27 05:58:56 MET DST 1998
1706	//------------------------------------------------------------
1707	// @function
1708
1709	//!@{
1710	void
1711	present(void)
1712	{
1713	cout << "##################################################\n"
1714	<< SIGNATURE << '\n' << '\n'
1715	<< "Processor of the reflector output\n"
1716	<< "J C Gonzalez, Jun 1998\n"
1717	<< "##################################################\n\n"
1718	<< flush ;
1719	}
1720	//!@}
1721
1722
1723	//!-----------------------------------------------------------
1724	// @name usage
1725	//
1726	// @desc show help
1727	//
1728	// @date Tue Dec 15 16:23:30 MET 1998
1729	//------------------------------------------------------------
1730	// @function
1731
1732	//!@{
1733	void
1734	usage(void)
1735	{
1736	present();
1737	cout << "\nusage ::\n\n"
1738	<< "\t camera "
1739	<< " [ -@ paramfile ] "
1740	<< " [ -h ] "
1741	<< "\n\n or \n\n"
1742	<< "\t camera < paramfile"
1743	<< "\n\n";
1744	exit(0);
1745	}
1746	//!@}
1747
1748
1749	//!-----------------------------------------------------------
1750	// @name log
1751	//
1752	// @desc function to send log information
1753	//
1754	// @var funct Name of the caller function
1755	// @var fmt Format to be used (message)
1756	// @var ... Other information to be shown
1757	//
1758	// @date Sat Jun 27 05:58:56 MET DST 1998
1759	//------------------------------------------------------------
1760	// @function
1761
1762	//!@{
1763	void
1764	log(const char funct, char fmt, ...)
1765	{
1766	va_list args;
1767
1768	// Display the name of the function that called error
1769	printf("[%s]: ", funct);
1770
1771	// Display the remainder of the message
1772	va_start(args, fmt);
1773	vprintf(fmt, args);
1774	va_end(args);
1775	}
1776	//!@}
1777
1778
1779	//!-----------------------------------------------------------
1780	// @name error
1781	//
1782	// @desc function to send an error message, and abort the program
1783	//
1784	// @var funct Name of the caller function
1785	// @var fmt Format to be used (message)
1786	// @var ... Other information to be shown
1787	//
1788	// @date Sat Jun 27 05:58:56 MET DST 1998
1789	//------------------------------------------------------------
1790	// @function
1791
1792	//!@{
1793	void
1794	error(const char funct, char fmt, ...)
1795	{
1796	va_list args;
1797
1798	// Display the name of the function that called error
1799	fprintf(stderr, "ERROR in %s: ", funct);
1800
1801	// Display the remainder of the message
1802	va_start(args, fmt);
1803	vfprintf(stderr, fmt, args);
1804	va_end(args);
1805
1806	perror(funct);
1807
1808	abort();
1809	}
1810	//!@}
1811
1812
1813	//!-----------------------------------------------------------
1814	// @name isA
1815	//
1816	// @desc returns TRUE(FALSE), if the flag is(is not) the given
1817	//
1818	// @var s1 String to be searched
1819	// @var flag Flag to compare with string s1
1820	// @return TRUE: both strings match; FALSE: oth.
1821	//
1822	// @date Wed Jul 8 15:25:39 MET DST 1998
1823	//------------------------------------------------------------
1824	// @function
1825
1826	//!@{
1827	int
1828	isA( char * s1, const char * flag ) {
1829	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
1830	}
1831	//!@}
1832
1833
1834	//!-----------------------------------------------------------
1835	// @name read_ct_file
1836	//
1837	// @desc read CT definition file
1838	//
1839	// @date Sat Jun 27 05:58:56 MET DST 1998
1840	//------------------------------------------------------------
1841	// @function
1842
1843	//!@{
1844	void
1845	read_ct_file(void)
1846	{
1847	char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
1848	char token[ITEM_MAX_LENGTH]; //@< a single token
1849	int i, j; //@< dummy counters
1850
1851	log( "read_ct_file", "start.\n" );
1852
1853	ifstream ctin ( ct_filename );
1854
1855	if ( ctin.bad() )
1856	error( "read_ct_file",
1857	"Cannot open CT def. file: %s\n", ct_filename );
1858
1859	// loop till the "end" directive is reached
1860
1861	while (!ctin.eof()) {
1862
1863	// get line from stdin
1864
1865	ctin.getline(line, LINE_MAX_LENGTH);
1866
1867	// look for each item at the beginning of the line
1868
1869	for (i=0; i<=define_mirrors; i++)
1870	if (strstr(line, CT_ITEM_NAMES[i]) == line)
1871	break;
1872
1873	// if it is not a valid line, just ignore it
1874
1875	if (i == define_mirrors+1)
1876	continue;
1877
1878	// case block for each directive
1879
1880	switch ( i ) {
1881
1882	case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
1883
1884	// get focal distance
1885
1886	sscanf(line, "%s %d", token, &ct_Type);
1887
1888	log( "read_ct_file", "<Type of Telescope>: %s\n",
1889	((ct_Type==0) ? "CT1" : "MAGIC") );
1890
1891	break;
1892
1893	case focal_distance: // <focal distance> [cm]
1894
1895	// get focal distance
1896
1897	sscanf(line, "%s %f", token, &ct_Focal_mean);
1898
1899	log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
1900
1901	break;
1902
1903	case focal_std: // s(focal distance) [cm]
1904
1905	// get focal distance
1906
1907	sscanf(line, "%s %f", token, &ct_Focal_std);
1908
1909	log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
1910
1911	break;
1912
1913	case point_spread: // <point spread> [cm]
1914
1915	// get point spread
1916
1917	sscanf(line, "%s %f", token, &ct_PSpread_mean);
1918
1919	log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
1920
1921	break;
1922
1923	case point_std: // s(point spread) [cm]
1924
1925	// get point spread
1926
1927	sscanf(line, "%s %f", token, &ct_PSpread_std);
1928
1929	log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
1930
1931	break;
1932
1933	case adjustment_dev: // s(adjustment_dev) [cm]
1934
1935	// get point spread
1936
1937	sscanf(line, "%s %f", token, &ct_Adjustment_std);
1938
1939	log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
1940
1941	break;
1942
1943	case black_spot: // radius of the black spot in the center [cm]
1944
1945	// get black spot radius
1946
1947	sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
1948
1949	log( "read_ct_file", "Radius of the black spots: %f cm\n",
1950	ct_BlackSpot_rad);
1951
1952	break;
1953
1954	case r_mirror: // radius of the mirrors [cm]
1955
1956	// get radius of mirror
1957
1958	sscanf(line, "%s %f", token, &ct_RMirror);
1959
1960	log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
1961
1962	break;
1963
1964	case n_mirrors: // number of mirrors
1965
1966	// get the name of the output_file from the line
1967
1968	sscanf(line, "%s %d", token, &ct_NMirrors);
1969
1970	log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
1971
1972	break;
1973
1974	case camera_width: // camera width [cm]
1975
1976	// get the name of the ct_file from the line
1977
1978	sscanf(line, "%s %f", token, &ct_CameraWidth);
1979
1980	log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
1981
1982	break;
1983
1984	case n_pixels: // number of pixels
1985
1986	// get the name of the output_file from the line
1987
1988	sscanf(line, "%s %d", token, &ct_NPixels);
1989
1990	log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
1991
1992	break;
1993
1994	case n_centralpixels: // number of central pixels
1995
1996	// get the name of the output_file from the line
1997
1998	sscanf(line, "%s %d", token, &ct_NCentralPixels);
1999
2000	log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
2001
2002	break;
2003
2004	case n_gappixels: // number of gap pixels
2005
2006	// get the name of the output_file from the line
2007
2008	sscanf(line, "%s %d", token, &ct_NGapPixels);
2009
2010	log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
2011
2012	break;
2013
2014	case pixel_width: // pixel width [cm]
2015
2016	// get the name of the ct_file from the line
2017
2018	sscanf(line, "%s %f", token, &ct_PixelWidth);
2019
2020	ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
2021	ct_PixelWidth_corner_2_corner_half =
2022	ct_PixelWidth_corner_2_corner * 0.50;
2023	ct_Apot = ct_PixelWidth / 2;
2024	ct_2Apot = ct_Apot * 2.0;
2025
2026	log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2027
2028	break;
2029
2030	case define_mirrors: // read table with the parameters of the mirrors
2031
2032	log( "read_ct_file", "Table of mirrors data:\n" );
2033
2034	// check whether the number of mirrors was already set
2035
2036	if ( ct_NMirrors == 0 )
2037	error( "read_ct_file", "NMirrors was not set.\n" );
2038
2039	// allocate memory for paths list
2040
2041	log( "read_ct_file", "Allocating memory for ct_data\n" );
2042
2043	ct_data = new float*[ct_NMirrors];
2044
2045	for (i=0; i<ct_NMirrors; i++)
2046	ct_data[i] = new float[CT_NDATA];
2047
2048	// read data
2049
2050	log( "read_ct_file", "Reading mirrors data...\n" );
2051
2052	for (i=0; i<ct_NMirrors; i++)
2053	for (j=0; j<CT_NDATA; j++)
2054	ctin >> ct_data[i][j];
2055
2056	break;
2057
2058	} // switch ( i )
2059
2060	} // end while
2061
2062	// end
2063
2064	log( "read_ct_file", "done.\n" );
2065
2066	return;
2067	}
2068	//!@}
2069
2070
2071	//!-----------------------------------------------------------
2072	// @name read_pixels
2073	//
2074	// @desc read pixels data
2075	//
2076	// @date Fri Mar 12 16:33:34 MET 1999
2077	//------------------------------------------------------------
2078	// @function
2079
2080	//!@{
2081	void
2082	read_pixels(struct camera *pcam)
2083	{
2084	ifstream qefile;
2085	char line[LINE_MAX_LENGTH];
2086	int n, i, j, k;
2087	float qe;
2088
2089	//------------------------------------------------------------
2090	// first, pixels' coordinates
2091
2092	pcam->inumpixels = ct_NPixels;
2093	pcam->inumcentralpixels = ct_NCentralPixels;
2094	pcam->inumgappixels = ct_NGapPixels;
2095	pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
2096	pcam->dpixdiameter_cm = ct_PixelWidth;
2097
2098	// initialize pixel numbers
2099
2100	for ( i=0; i<PIX_ARRAY_SIDE; ++i )
2101	for ( j=0; j<PIX_ARRAY_SIDE; ++j )
2102	pixels[i][j][PIXNUM] = -1;
2103
2104	pixary = new float* [2*ct_NCentralPixels];
2105	for ( i=0; i<2*ct_NCentralPixels; ++i )
2106	pixary[i] = new float[2];
2107
2108	pixneig = new int* [ct_NCentralPixels];
2109	for ( i=0; i<ct_NCentralPixels; ++i ) {
2110	pixneig[i] = new int[6];
2111	for ( j=0; j<6; ++j )
2112	pixneig[i][j] = -1;
2113	}
2114
2115	npixneig = new int[ct_NCentralPixels];
2116	for ( i=0; i<ct_NCentralPixels; ++i )
2117	npixneig[i] = 0;
2118
2119	// generate all coordinates
2120
2121	igen_pixel_coordinates(pcam);
2122
2123	// transfer coordinates to the working arrays for
2124	// the central pixels
2125
2126	for(k=0; k<ct_NCentralPixels; k++){
2127
2128	i = (int) pcam->di[k];
2129	j = (int) pcam->dj[k];
2130
2131	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k;
2132	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k];
2133	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k];
2134
2135	pixary[k][0] = pcam->dxc[k];
2136	pixary[k][1] = pcam->dyc[k];
2137	}
2138
2139	// calculate tables of neighbours
2140
2141	#ifdef __DEBUG__
2142	for ( n=0 ; n<ct_NPixels ; ++n ) {
2143	cout << "Para el pixel " << n << ": ";
2144	for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2145	if ( pixels_are_neig(n,i) == TRUE ) {
2146	pixneig[n][npixneig[n]] = i;
2147	pixneig[i][npixneig[i]] = n;
2148	cout << i << ' ';
2149	++npixneig[n];
2150	++npixneig[i];
2151	}
2152	}
2153	cout << endl << flush;
2154	}
2155	#else // ! __DEBUG__
2156	for ( n=0 ; n<ct_NCentralPixels ; ++n )
2157	for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2158	if ( pixels_are_neig(n,i) == TRUE ) {
2159	pixneig[n][npixneig[n]] = i;
2160	pixneig[i][npixneig[i]] = n;
2161	++npixneig[n];
2162	++npixneig[i];
2163	}
2164	#endif // ! __DEBUG__
2165
2166	#ifdef __DEBUG__
2167	for ( n=0 ; n<ct_NPixels ; ++n ) {
2168	cout << n << ':';
2169	for ( j=0; j<npixneig[n]; ++j)
2170	cout << ' ' << pixneig[n][j];
2171	cout << endl << flush;
2172	}
2173	#endif // __DEBUG__
2174
2175	//------------------------------------------------------------
2176	// second, pixels' QE
2177
2178	// try to open the file
2179
2180	log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE);
2181
2182	qefile.open( QE_FILE );
2183
2184	// if it is wrong or does not exist, exit
2185
2186	if ( qefile.bad() )
2187	error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2188
2189	// read file
2190
2191	log("read_pixels", "Reading data . . .\n");
2192
2193	i=-1;
2194
2195	while ( ! qefile.eof() ) {
2196
2197	// get line from the file
2198
2199	qefile.getline(line, LINE_MAX_LENGTH);
2200
2201	// skip if comment
2202
2203	if ( *line == '#' )
2204	continue;
2205
2206	// if it is the first valid value, it is the number of QE data points
2207
2208	if ( i < 0 ) {
2209
2210	// get the number of datapoints
2211
2212	sscanf(line, "%d", &pointsQE);
2213
2214	// allocate memory for the table of QEs
2215
2216	QE = new float ** [ct_NPixels];
2217
2218	for ( i=0; i<ct_NPixels; ++i ) {
2219	QE[i] = new float * [2];
2220	QE[i][0] = new float[pointsQE];
2221	QE[i][1] = new float[pointsQE];
2222	}
2223
2224	QElambda = new float [pointsQE];
2225
2226	for ( i=0; i<pointsQE; ++i ) {
2227	qefile.getline(line, LINE_MAX_LENGTH);
2228	sscanf(line, "%f", &QElambda[i]);
2229	}
2230
2231	i=0;
2232
2233	continue;
2234	}
2235
2236	// get the values (num-pixel, num-datapoint, QE-value)
2237
2238	sscanf(line, "%d %d %f", &i, &j, &qe);
2239
2240	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2241	((j-1) < pointsQE) && ((j-1) > -1) ) {
2242	QE[i-1][0][j-1] = QElambda[j-1];
2243	QE[i-1][1][j-1] = qe;
2244	}
2245
2246	}
2247
2248	// close file
2249
2250	qefile.close();
2251
2252	// end
2253
2254	log("read_pixels", "Done.\n");
2255
2256	}
2257	//!@}
2258
2259
2260	//!-----------------------------------------------------------
2261	// @name pixels_are_neig
2262	//
2263	// @desc check whether two pixels are neighbours
2264	//
2265	// @var pix1 Number of the first pixel
2266	// @var pix2 Number of the second pixel
2267	// @return TRUE: both pixels are neighbours; FALSE: oth.
2268	//
2269	// @date Wed Sep 9 17:58:37 MET DST 1998
2270	//------------------------------------------------------------
2271	// @function
2272
2273	//!@{
2274	int
2275	pixels_are_neig(int pix1, int pix2)
2276	{
2277	if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2278	SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2279	> ct_PixelWidth_corner_2_corner )
2280	return ( FALSE );
2281	else
2282	return ( TRUE );
2283	}
2284	//!@}
2285
2286	//!-----------------------------------------------------------
2287	// @name igen_pixel_coordinates
2288	//
2289	// @desc generate the pixel center coordinates
2290	//
2291	// @var *pcam structure camera containing all the
2292	// camera information
2293	// @return total number of pixels
2294	//
2295	// DP
2296	//
2297	// @date Thu Oct 14 10:41:03 CEST 1999
2298	//------------------------------------------------------------
2299	// @function
2300
2301	//!@{
2302	/****** igen_pixel_coordinates() *******************************/
2303
2304	int igen_pixel_coordinates(struct camera *pcam) {
2305	/* generate pixel coordinates, return value is number of pixels */
2306
2307	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2308	float fsegment_fract;
2309	double dtsize;
2310	double dhsize;
2311	double dpsize;
2312	double dxfirst_pix;
2313	double dyfirst_pix;
2314	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
2315	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
2316
2317
2318	double dstartx, dstarty; /* for the gap pixels and outer pixels */
2319	int j, nrow;
2320
2321	dpsize = pcam->dpixdiameter_cm;
2322	dtsize = dpsize * sqrt(3.) / 2.;
2323	dhsize = dpsize / 2.;
2324
2325	/* Loop over central pixels to generate co-ordinates */
2326
2327	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2328
2329	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2330
2331	pcam->dpixsizefactor[ipixno] = 1.;
2332
2333	in = 0;
2334
2335	i = 0;
2336	itot_inside_ring = 0;
2337	iring_no = 0;
2338
2339	/* Calculate the number of pixels out to and including the ring containing pixel number */
2340	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2341
2342	while (itot_inside_ring == 0){
2343
2344	iN = 3(i(i+1)) + 1;
2345
2346	if (ipixno <= iN){
2347	iring_no = i;
2348	itot_inside_ring = iN;
2349	}
2350
2351	i++;
2352	}
2353
2354
2355	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2356
2357	ipix_in_ring = 0;
2358	for (i = 0; i < iring_no; ++i){
2359
2360	ipix_in_ring = ipix_in_ring + 6;
2361	}
2362
2363	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2364	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2365	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2366
2367	isegment = 0;
2368	fsegment_fract = 0.;
2369	if (iring_no > 0) {
2370
2371	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2372
2373	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2374
2375	}
2376
2377	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2378	/* pixel width (flat to flat) dpsize. */
2379
2380	dxfirst_pix = dpsize*iring_no;
2381	dyfirst_pix = 0.;
2382
2383	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2384
2385	ddxseg1 = - dhsize*iring_no;
2386	ddyseg1 = dtsize*iring_no;
2387	ddxseg2 = -dpsize*iring_no;
2388	ddyseg2 = 0.;
2389	ddxseg3 = ddxseg1;
2390	ddyseg3 = -ddyseg1;
2391	ddxseg4 = -ddxseg1;
2392	ddyseg4 = -ddyseg1;
2393	ddxseg5 = -ddxseg2;
2394	ddyseg5 = 0.;
2395	ddxseg6 = -ddxseg1;
2396	ddyseg6 = ddyseg1;
2397
2398	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2399	/* anti-clockwise around the ring by adding the segment to segment vectors. */
2400
2401	switch (isegment) {
2402
2403	case 0:
2404
2405	pcam->dxc[ipixno-1] = 0.;
2406	pcam->dyc[ipixno-1] = 0.;
2407
2408	case 1:
2409	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2410	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2411
2412	break;
2413
2414	case 2:
2415
2416	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2417	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2418
2419	break;
2420
2421	case 3:
2422
2423	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
2424	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
2425
2426	break;
2427
2428	case 4:
2429
2430	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
2431	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
2432
2433	break;
2434
2435	case 5:
2436
2437	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
2438	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2439
2440	break;
2441
2442	case 6:
2443
2444	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2445	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2446
2447	break;
2448
2449	default:
2450
2451	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2452	return(0);
2453
2454	} /* end switch */
2455
2456	} /* end for */
2457
2458	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2459	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2460
2461	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2462
2463	j = pcam->inumcentralpixels;
2464
2465	for(i=0; i<pcam->inumgappixels; i=i+6){
2466	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2467	pcam->dyc[j + i ] = dstarty;
2468	pcam->dpixsizefactor[j + i] = 1.;
2469	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2470	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2471	pcam->dpixsizefactor[j + i + 1] = 1.;
2472	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2473	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2474	pcam->dpixsizefactor[j + i+ 2] = 1.;
2475	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2476	pcam->dyc[j + i + 3] = dstarty;
2477	pcam->dpixsizefactor[j + i+ 3] = 1.;
2478	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2479	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2480	pcam->dpixsizefactor[j + i+ 4] = 1.;
2481	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2482	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2483	pcam->dpixsizefactor[j + i + 5] = 1.;
2484	} /* end for */
2485	} /* end if */
2486
2487	/* generate positions of the outer pixels */
2488
2489	if( pcam->inumbigpixels > 0 ){
2490
2491	j = pcam->inumcentralpixels + pcam->inumgappixels;
2492
2493	for(i=0; i<pcam->inumbigpixels; i++){
2494	pcam->dpixsizefactor[j + i] = 2.;
2495	}
2496
2497	in = 0;
2498
2499	nrow = (int) ceil(dstartx / 2. / dpsize);
2500
2501	while(in < pcam->inumbigpixels){
2502
2503	pcam->dxc[j + in] = dstartx + dpsize;
2504	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2505	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2506	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2507	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2508	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2509	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2510	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2511	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2512	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2513	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2514	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2515	for(i=1; i<nrow; i++){
2516	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2517	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2518	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2519	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2520	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2521	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2522	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2523	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2524	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
2525	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
2526	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
2527	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
2528	}
2529	in = in + 6 * nrow;
2530	dstartx = dstartx + 2. * dpsize;
2531	nrow = nrow + 1;
2532	} /* end while */
2533
2534	} /* end if */
2535
2536	/* generate the ij coordinates */
2537
2538	for(i=0; i<pcam->inumpixels; i++){
2539	pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize;
2540	pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60;
2541
2542	// fprintf(stdout, "%d %f %f %f %f %f\n",
2543	// i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i],
2544	// pcam->dpixsizefactor[i]);
2545
2546	}
2547
2548	return(pcam->inumpixels);
2549
2550	}
2551	//!@}
2552
2553	//!-----------------------------------------------------------
2554	// @name bpoint_is_in_pix
2555	//
2556	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
2557	//
2558	// @var *pcam structure camera containing all the
2559	// camera information
2560	// @var dx, dy point coordinates in centimeters
2561	// @var ipixnum pixel number (starting at 0)
2562	// @return TRUE if the point is inside the pixel, FALSE otherwise
2563	//
2564	// DP
2565	//
2566	// @date Thu Oct 14 16:59:04 CEST 1999
2567	//------------------------------------------------------------
2568	// @function
2569
2570	//!@{
2571
2572	/****** bpoint_is_in_pix() *************************************/
2573
2574	int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
2575	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
2576	/* the pixel is assumed to be a "closed set" */
2577
2578	double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
2579	double c, xx, yy; /* auxiliary variable */
2580
2581	b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
2582	a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
2583	c = 1. - 1./sqrt(3.);
2584	if((ipixnum < 0)\|\|(ipixnum >= pcam->inumpixels)){
2585	fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
2586	fprintf(stderr, "Exiting.\n");
2587	exit(203);
2588	}
2589	xx = dx - pcam->dxc[ipixnum];
2590	yy = dy - pcam->dyc[ipixnum];
2591
2592	if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) \|\|
2593	((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
2594	return(TRUE); /* point is inside */
2595	}
2596	else{
2597	return(FALSE); /* point is outside */
2598	}
2599	}
2600
2601	//!@}
2602
2603	//------------------------------------------------------------
2604	// @name dist_r_P
2605	//
2606	// @desc distance straight line r - point P
2607	//
2608	// @date Sat Jun 27 05:58:56 MET DST 1998
2609	// @function @code
2610	//------------------------------------------------------------
2611	// dist_r_P
2612	//
2613	// distance straight line r - point P
2614	//------------------------------------------------------------
2615
2616	float
2617	dist_r_P(float a, float b, float c,
2618	float l, float m, float n,
2619	float x, float y, float z)
2620	{
2621	return (
2622	sqrt((SQR((a-x)m-(b-y)l) +
2623	SQR((b-y)n-(c-z)m) +
2624	SQR((c-z)l-(a-x)n))/
2625	(SQR(l)+SQR(m)+SQR(n))
2626	)
2627	);
2628	}
2629	// @endcode
2630
2631
2632	//=------------------------------------------------------------
2633	//!@subsection Log of this file.
2634
2635	//!@{
2636	//
2637	// $Log: not supported by cvs2svn $
2638	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
2639	// This the starting point for CVS controlled further developments of the
2640	// camera program. The program was originally written by Jose Carlos.
2641	// But here you can find a "rootified" version to the program. This means
2642	// that there is no hbook stuff in it now. Also the output of the
2643	// program changed to the MagicRawDataFormat.
2644	//
2645	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
2646	//
2647	// In the following you can see the README file of that version:
2648	//
2649	// ==================================================
2650	//
2651	// Fri Oct 22 1999 D.P.
2652	//
2653	// The MAGIC Monte Carlo System
2654	//
2655	// Camera Simulation Programme
2656	// ---------------------------
2657	//
2658	// 1) Description
2659	//
2660	// This version is the result of the fusion of H.K.'s
2661	// root_camera which is described below (section 2)
2662	// and another version by D.P. which had a few additional
2663	// useful features.
2664	//
2665	// The version compiles under Linux with ROOT 2.22 installed
2666	// (variable ROOTSYS has to be set).
2667	//
2668	// Compile as before simply using "make" in the root_camera
2669	// directory.
2670	//
2671	// All features of H.K.'s root_camera were retained.
2672	//
2673	// Additional features of this version are:
2674	//
2675	// a) HBOOK is no longer used and all references are removed.
2676	//
2677	// b) Instead of HBOOK, the user is given now the possibility of
2678	// having Diagnostic data in ROOT format as a complement
2679	// to the ROOT Raw data.
2680	//
2681	// This data is written to the file which is determined by
2682	// the new input parameter "diag_file" in the camera parameter
2683	// file.
2684	//
2685	// All source code file belonging to this part have filenames
2686	// starting with "MDiag".
2687	//
2688	// The user can read the output file using the following commands
2689	// in an interactive ROOT session:
2690	//
2691	// root [0] .L MDiag.so
2692	// root [1] new TFile("diag.root");
2693	// root [2] new TTreeViewer("T");
2694	//
2695	// This brings up a viewer from which all variables of the
2696	// TTree can be accessed and histogrammed. This example
2697	// assumes that you have named the file "diag.root", that
2698	// you are using ROOT version 2.22 or later and that you have
2699	// the shared object library "MDiag.so" which is produced
2700	// by the Makefile along with the executable "camera".
2701	//
2702	// ! The contents of the so-called diag file is not yet fixed.
2703	// ! At the moment it is what J.C.G. used to put into the HBOOK
2704	// ! ntuple. In future versions the moments calculation can be
2705	// ! removed and the parameter list be modified correspondingly.
2706	//
2707	// c) Now concatenated reflector files can be read. This is useful
2708	// if you have run the reflector with different parameters but
2709	// you want to continue the analysis with all reflector data
2710	// going into ONE ROOT outputfile.
2711	//
2712	// The previous camera version contained a bug which made reading
2713	// of two or more concatenated reflector files impossible.
2714	//
2715	// d) The reflector output format was changed. It is now version
2716	// 0.4 .
2717	// The change solely consists in a shortening of the flag
2718	// definition in the file
2719	//
2720	// include-MC/MCCphoton.hxx
2721	//
2722	// ! IF YOU WANT TO READ REFLECTOR FORMAT 0.3, you can easily
2723	// ! do so by recompiling camera with the previous version of
2724	// ! include-MC/MCCphoton.hxx.
2725	//
2726	// The change was necessary for saving space and better
2727	// debugging. From now on, this format can be frozen.
2728	//
2729	// ! For producing reflector output in the new format, you
2730	// ! of course have to recompile your reflector with the
2731	// ! new include-MC/MCCphoton.hxx .
2732	//
2733	// e) A first version of the pixelization with the larger
2734	// outer pixels is implemented. THIS IS NOT YET FULLY
2735	// TESTED, but first rough tests show that it works
2736	// at least to a good approximation.
2737	//
2738	// The present version implements the camera outline
2739	// with 18 "gap-pixels" and 595 pixels in total as
2740	// shown in
2741	//
2742	// http://sarastro.ifae.es/internal/home/hardware/camera/numbering.ps
2743	//
2744	// This change involved
2745	//
2746	// (i) The file pixels.dat is no longer needed. Instead
2747	// the coordinates are generated by the program itself
2748	// (takes maybe 1 second). In the file
2749	//
2750	// pixel-coords.txt
2751	//
2752	// in the same directory as this README, you find a list
2753	// of the coordinates generated by this new routine. It
2754	// has the format
2755	//
2756	// number i j x y size-factor
2757	//
2758	// where i and j are J.C.G.'s so called biaxis hexagonal
2759	// coordinates (for internal use) and x and y are the
2760	// coordinates of the pixel centers in the standard camera
2761	// coordinate system in units of centimeters. The value
2762	// of "size-factor" determines the linear size of the pixel
2763	// relative to the central pixels.
2764	//
2765	// (ii) The magic.def file has two additional parameters
2766	// which give the number of central pixels and the
2767	// number of gap pixels
2768	//
2769	// (iii) In camera.h and camera.cxx several changes were
2770	// necessary, among them the introduction of several
2771	// new functions
2772	//
2773	// The newly suggested outline with asymmetric Winston cones
2774	// will be implemented in a later version.
2775	//
2776	// f) phe files can no longer be read since this contradicts
2777	// our philosophy that the analysis should be done with other
2778	// programs like e.g. EVITA and not with "camera" itself.
2779	// This possibility was removed.
2780	//
2781	// g) ROOT is no longer invoked with an interactive interface.
2782	// In this way, camera can better be run as a batch program and
2783	// it uses less memory.
2784	//
2785	// h) small changes concerning the variable "t_chan" were necessary in
2786	// order to avoid segmentation faults: The variable is used as an
2787	// index and it went sometimes outside the limits when camera
2788	// was reading proton data. This is because the reflector files
2789	// don't contain the photons in a chronological order and also
2790	// the timespread can be considerably longer that the foreseen
2791	// digitisation timespan. Please see the source code of camera.cxx
2792	// round about line 1090.
2793	//
2794	// j) several unused variables were removed, a few warning messages
2795	// occur when you compile camera.cxx but these can be ignored at
2796	// the moment.
2797	//
2798	// In general the program is of course not finished. It still needs
2799	// debugging, proper trigger simulation, simulation of the asymmetric
2800	// version of the outer pixels, proper NSB simulation, adaption of
2801	// the diag "ntuple" contents to our need and others small improvements.
2802	//
2803	// In the directory rfl-files there is now a file in reflector format 0.4
2804	// containing a single event produced by the starfiled adder. It has
2805	// a duration of 30 ns and represents the region around the Crab Nebula.
2806	// Using the enclosed input parameter file, camera should process this
2807	// file without problems.
2808	//
2809	// 2) The README for the previous version of root_camera
2810	//
2811	// README for a preliminary version of the
2812	// root_camera program.
2813	//
2814	// root_camera is based on the program "camera"of Jose Carlos
2815	// Gonzalez. It was changed in the way that only the pixelisation
2816	// and the distibution of the phe to the FADCs works in a
2817	// first version.
2818	//
2819	// Using the #undef command most possibilities of the orignal
2820	// program are switched of.
2821	//
2822	// The new parts are signed by
2823	//
2824	// - ROOT or __ROOT__
2825	// nearly all important codelines for ROOT output are enclosed
2826	// in structures like
2827	// #ifdef __ROOT__
2828	//
2829	// code
2830	//
2831	// #endif // __ROOT__
2832	//
2833	// In same case the new lines are signed by a comment with the word
2834	// ROOT in it.
2835	//
2836	// For timing of the pulse some variable names are changed.
2837	// (t0, t1, t --> t_ini, t_fin, t_1st, t_chan,...)
2838	// Look also for this changes.
2839	//
2840	// For the new root-file is also a change in readparm-files
2841	//
2842	//
2843	// - __DETAIL_TRIGGER__
2844	//
2845	// This is for the implementation of the current work on trigger
2846	// studies. Because the class MTrigger is not well documented it
2847	// isn´t a part of this tar file. Only a dummy File exists.
2848	//
2849	//
2850	//
2851	// With all files in the archive, the root_camera program should run.
2852	//
2853	// A reflector file is in the directory rfl-files
2854	//
2855	// ==================================================
2856	//
2857	// From now on, use CVS for development!!!!
2858	//
2859	//
2860	//
2861	// Revision 1.3 1999/10/22 15:01:28 petry
2862	// version sent to H.K. and N.M. on Fri Oct 22 1999
2863	//
2864	// Revision 1.2 1999/10/22 09:44:23 petry
2865	// first synthesized version which compiles and runs without crashing;
2866	//
2867	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
2868	// first synthesised version
2869	//
2870	// Revision 1.13 1999/03/15 14:59:05 gonzalez
2871	// camera-1_1
2872	//
2873	// Revision 1.12 1999/03/02 09:56:10 gonzalez
2874	// * empty log message *
2875	//
2876	//
2877	//!@}
2878
2879	//=EOF

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

Download in other formats: