Context Navigation

source: trunk/MagicSoft/Simulation/Detector/TimeCam/timecam.cxx@ 372

Visit:

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

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

Download in other formats: