Context Navigation

timecam.cxx@ 469

Visit:

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

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

Context Navigation

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

Download in other formats: