Context Navigation

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

Visit:

Last change on this file since 357 was 357, checked in by harald, 26 years ago
small changes to run the timecam program.
File size: 74.4 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.2 $
24	// $Author: harald $
25	// $Date: 2000-02-09 12:45:28 $
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	if ( wl < 1.0 )
1171	break;
1172
1173	if ( (wl > 600.0) \|\| (wl < 290.0) )
1174	break;
1175
1176	// check if photon is inside outermost camera radius
1177
1178	if(sqrt(cxcx + cycy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)){
1179
1180	// read next CPhoton
1181	if ( Data_From_STDIN )
1182	cin.read( flag, SIZE_OF_FLAGS );
1183	else
1184	inputfile.read ( flag, SIZE_OF_FLAGS );
1185
1186	// go to beginning of loop, the photon is lost
1187	continue;
1188
1189	}
1190
1191	// cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
1192
1193	ci = floor( (cx - cy*COS60/SIN60)/ ct_2Apot + 0.5);
1194	cj = floor( (cy/SIN60) / ct_2Apot + 0.5);
1195
1196	ici = (int)(ci);
1197	icj = (int)(cj);
1198
1199	iici = ici+PIX_ARRAY_HALF_SIDE;
1200	iicj = icj+PIX_ARRAY_HALF_SIDE;
1201
1202	// is it inside the array?
1203
1204	if ( (iici > 0) && (iici < PIX_ARRAY_SIDE) &&
1205	(iicj > 0) && (iicj < PIX_ARRAY_SIDE) ) {
1206
1207	// try to put into pixel
1208
1209	// obtain the pixel number for this photon
1210
1211	nPMT = (int)
1212	pixels[ici+PIX_ARRAY_HALF_SIDE][icj+PIX_ARRAY_HALF_SIDE][PIXNUM];
1213
1214	}
1215	else{
1216
1217	nPMT = -1;
1218
1219	}
1220
1221	// check if outside the central camera
1222
1223	if ( (nPMT < 0) \|\| (nPMT >= ct_NCentralPixels) ) {
1224
1225	// check the outer pixels
1226	nPMT = -1;
1227
1228	for(i=ct_NCentralPixels; i<ct_NPixels; i++){
1229	if( bpoint_is_in_pix( cx, cy, i, &cam) ){
1230	nPMT = i;
1231	break;
1232	}
1233	}
1234
1235	if(nPMT==-1){// the photon is in none of the pixels
1236
1237	// read next CPhoton
1238	if ( Data_From_STDIN )
1239	cin.read( flag, SIZE_OF_FLAGS );
1240	else
1241	inputfile.read ( flag, SIZE_OF_FLAGS );
1242
1243	// go to beginning of loop, the photon is lost
1244	continue;
1245	}
1246
1247	}
1248
1249	#ifdef __QE__
1250
1251	//!@' @#### QE simulation.
1252	//@'
1253
1254	//+++
1255	// QE simulation
1256	//---
1257
1258	// find data point to be used in Lagrange interpolation (-> k)
1259
1260	qeptr = (float **)QE[nPMT];
1261
1262	FindLagrange(qeptr,k,wl);
1263
1264	// if random > quantum efficiency, reject it
1265
1266	qe = Lagrange(qeptr,k,wl) / 100.0;
1267
1268	// fprintf(stdout, "%f\n", qe);
1269
1270	if ( RandomNumber > qe ) {
1271
1272	// read next CPhoton
1273	if ( Data_From_STDIN )
1274	cin.read( flag, SIZE_OF_FLAGS );
1275	else
1276	inputfile.read ( flag, SIZE_OF_FLAGS );
1277
1278	// go to beginning of loop
1279	continue;
1280
1281	}
1282
1283	#endif // __QE__
1284
1285	//+++
1286	// Cphoton is accepted
1287	//---
1288
1289	ncph_out++ ;
1290
1291	// increase the number of Cphs. in the PMT, i.e.,
1292	// increase in one unit the counter of the photons
1293	// stored in the pixel nPMT
1294
1295	fnpix[nPMT] += 1.0;
1296
1297	#ifdef __DETAIL_TRIGGER__
1298	//
1299	// fill the Trigger class with this phe
1300	//
1301	//
1302	Trigger.Fill( nPMT, ( t - t_first ) ) ;
1303
1304	// fadc.Fill( nPMT, ( t - t_first ), Trigger.Fill( nPMT, ( t - t_first ) ) ) ;
1305
1306	#endif // __DETAIL_TRIGGER__
1307
1308	// read next CPhoton
1309
1310	if ( Data_From_STDIN )
1311	cin.read( flag, SIZE_OF_FLAGS );
1312	else
1313	inputfile.read ( flag, SIZE_OF_FLAGS );
1314
1315	} // end while, i.e. found end of event
1316
1317	if ( fmod ( nshow, 1000. ) == 0. )
1318	log(SIGNATURE,
1319	"End of this event: in: %d cphs(+%d). out: %d cphs(+%d). .\n",
1320	ncph_in, ntcph_in,
1321	ncph_out, ntcph_out);
1322
1323	// show number of photons
1324
1325	//cout << ncph_in << " photons read . . . " << endl << flush;
1326
1327	// skip it ?
1328
1329	for ( i=0; i<nSkip; ++i ) {
1330	if (Skip[i] == (nshow+ntshow)) {
1331	i = -1;
1332	break;
1333	}
1334	}
1335
1336	// if after the previous loop, the exit value of i is -1
1337	// then the shower number is in the list of showers to be
1338	// skipped
1339
1340	if (i == -1) {
1341	log(SIGNATURE, "\t\tskipped!\n");
1342	continue;
1343	}
1344
1345	/*!@'
1346
1347	After reading all the Cherenkov photons for a given event,
1348	we have in the table of number of photons for each pixel
1349	only the 'raw' amount of Cherenkov photons @$n_p@$. Now, we
1350	should take this number as the mean value of the
1351	distribution of photons in that pixel @$p@$, following a
1352	Poisson distribution.
1353
1354	@[ n_p \equiv \mu_p @]
1355
1356	and with this number the amount of light coming from the
1357	shower is calculated @$\hat{n}_p@$.
1358
1359	Then, we calculate the amount of Night Sky Background we
1360	must introduce in that pixel @$p@$. We calculate this using
1361	again a Poisson distribution with mean @$\mu_\mathrm{NSB}@$
1362	(defined in the \|camera.h\| file). The value of
1363	@$\mu_\mathrm{NSB}@$ is obtained from measurements. With this
1364	value, the amount of photons @$\hat{n}_\mathrm{NSB}@$ coming
1365	from the Night Sky Background is calculated.
1366
1367	Finally, the amount of photons for that pixels is:
1368	@[ \hat{n}_p^\mathrm{final} = \hat{n}_p + \hat{n}_\mathrm{NSB} @]
1369
1370	*/
1371
1372	// after reading all the photons, our camera is filled
1373
1374	if ( Select_Energy ) {
1375	if (( mcevth.get_energy() < Select_Energy_le ) \|\|
1376	( mcevth.get_energy() > Select_Energy_ue )) {
1377	log(SIGNATURE, "select_energy: shower rejected.\n");
1378	continue;
1379	}
1380	}
1381
1382	#ifdef __NSB__
1383
1384	//!@' @#### NSB (Night Sky Background) simulation.
1385	//@'
1386
1387	//+++
1388	// NSB simulation
1389	//---
1390
1391	// add NSB "noise"
1392	// TO DO: make meanNSB an array and read the contents from a file!
1393
1394	// if ( simulateNSB )
1395	// for ( i=0; i<ct_NPixels; ++i )
1396	// fnpix[i] += (float)ignpoi( meanNSB );
1397	// old version of Jose Carlos
1398
1399	if ( simulateNSB) {
1400	//
1401	// loop over all pixels and scramble the number
1402	// of NSB photons to put in it. For the number use
1403	// a poison distribution with a mean calculated from meanNSB
1404	// and the TOTAL_TRIGGER_TIME
1405	//
1406	for ( Int_t nsbPix = 0 ; nsbPix < CAMERA_PIXELS ; nsbPix++ ) {
1407	//
1408	//
1409	meanPois = meanNSB * TOTAL_TRIGGER_TIME ;
1410
1411	// loop over the scrambled number if Photons in this pixels
1412
1413	for ( Int_t photNSB=0; photNSB<GenNSB.Poisson(meanPois);
1414	photNSB++){
1415	//
1416	// now scramble the time at that the photo electron of the
1417	// NSB photon is leaving the photo cathod
1418	//
1419
1420	timeNSB = GenNSB.Rndm() * TOTAL_TRIGGER_TIME ;
1421
1422	Trigger.FillNSB ( nsbPix, timeNSB ) ;
1423
1424	}
1425
1426
1427	}
1428
1429	}
1430	#endif // __NSB__
1431
1432	// if we should apply any kind of correction, do it here.
1433
1434	for ( i=0; i<ct_NPixels; ++i )
1435	fnpix[i] *= fCorrection;
1436
1437	#ifdef __DETAIL_TRIGGER__
1438	//
1439	// now the noise of the electronic
1440	// (preamps, optical transmission,..) is introduced.
1441	// This is done inside the class MTrigger by the method ElecNoise.
1442	//
1443	Trigger.ElecNoise() ;
1444
1445	//
1446	// look if in all the signals in the trigger signal branch
1447	// is a possible Trigger. Therefore we habe to diskriminate all
1448	// the simulated analog signals (Method Diskriminate in class
1449	// MTrigger). We look simultanously for the moments at which
1450	// there are more than TRIGGER_MULTI pixels above the
1451	// CHANNEL_THRESHOLD.
1452	//
1453
1454	McTrig->SetZeroLevel( Lev0 = (Short_t) Trigger.Diskriminate() ) ;
1455	Lev1 = Lev2 = 0 ;
1456
1457	//
1458	// Start the First Level Trigger simulation
1459	//
1460
1461	if ( Lev0 > 0 ) {
1462	McTrig->SetFirstLevel ( Lev1 = Trigger.FirstLevel() ) ;
1463	}
1464
1465	#endif // __DETAIL_TRIGGER__
1466
1467	#ifdef __ROOT__
1468
1469	//
1470	// Fill the header of this event
1471	//
1472
1473	Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ;
1474
1475	//
1476	// fill the MMcEvt with all information
1477	//
1478
1479	McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1480	mcevth.get_energy(),
1481	mcevth.get_theta(),
1482	mcevth.get_phi(),
1483	mcevth.get_core(),
1484	mcevth.get_coreX(),
1485	mcevth.get_coreY(),
1486	impactD,
1487	ulli, ulli,
1488	(UShort_t) ncph_in,
1489	ulli,
1490	(UShort_t) ncph_out ) ;
1491
1492	//
1493	// write it out to the file outfile
1494	//
1495
1496	EvtTree.Fill() ;
1497
1498	#endif // __ROOT__
1499
1500	//
1501	// if a first level trigger occurred, then
1502	// 1. do some other stuff (not implemented)
1503	// 2. start the gui tool
1504
1505	if ( Lev1 > 0 ) {
1506
1507	// fadc.Scan( Trigger.GetFirstLevelTime(0) ) ;
1508
1509	#ifdef __INTERAKTIV__
1510	Trigger.ShowSignal(McEvt) ;
1511	#endif
1512	}
1513
1514
1515
1516	#ifdef __ROOT__
1517	// clear all
1518	Evt->Clear() ;
1519	McEvt->Clear() ;
1520	McTrig->Clear() ;
1521	#endif // __ROOT__
1522
1523
1524	//++++++++++++++++++++++++++++++++++++++++++++++++++
1525	// at this point we have a camera full of
1526	// ph.e.s
1527	// we should first apply the trigger condition,
1528	// and if there's trigger, then clean the image,
1529	// calculate the islands statistics and the
1530	// other parameters of the image (Hillas' parameters
1531	// and so on).
1532	//--------------------------------------------------
1533
1534	#ifdef __DEBUG__
1535	printf("\n");
1536
1537	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1538
1539	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1540
1541	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1542
1543	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1544
1545	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1546	(int)pixels[ici][icj][PIXNUM],
1547	pixels[ici][icj][PIXX],
1548	pixels[ici][icj][PIXY],
1549	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1550
1551	}
1552	}
1553	}
1554	}
1555
1556	for (i=0; i<ct_NPixels; ++i) {
1557	printf("%d (%d): ", i, npixneig[i]);
1558	for (j=0; j<npixneig[i]; ++i)
1559	printf(" %d", pixneig[i][j]);
1560	printf("\n");
1561	}
1562
1563	#endif // __DEBUG__
1564
1565
1566	//!@' @#### Save data.
1567	//@'
1568
1569	//++++++++++++++++++++++++++++++++++++++++++++++++++
1570	// we now have all information we want
1571	// the only thing we must do now is writing it to
1572	// the output file
1573	//--------------------------------------------------
1574
1575	//++
1576	// save the image to the file
1577	//--
1578
1579	if ( Data_From_STDIN )
1580	cin.read( flag, SIZE_OF_FLAGS );
1581	else
1582	inputfile.read ( flag, SIZE_OF_FLAGS );
1583
1584	} // end while there is a next event
1585
1586	if( !isA( flag, FLAG_END_OF_RUN )){
1587	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1588	}
1589	else { // found end of run
1590	ntshow += nshow;
1591	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1592
1593	// if ( Data_From_STDIN )
1594	// cin.read( flag, SIZE_OF_FLAGS );
1595	// else
1596	// inputfile.read ( flag, SIZE_OF_FLAGS );
1597
1598	// huschel start here
1599
1600	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1601	log(SIGNATURE, "End of file . . .\n");
1602	still_in_loop = FALSE;
1603
1604	if ((! Data_From_STDIN) && (! inputfile.eof())){
1605
1606	// we have concatenated input files.
1607	// get signature of the next part and check it.
1608	// NOTE: this part repeats further up in the code;
1609	// if you change something here you probably want to change it
1610	// there as well
1611
1612	strcpy(Signature, REFL_SIGNATURE);
1613
1614	strcpy(sign, Signature);
1615
1616	inputfile.read( (char *)sign, strlen(Signature));
1617
1618	if (strcmp(sign, Signature) != 0) {
1619	cerr << "ERROR: Signature of .rfl file is not correct\n";
1620	cerr << '"' << sign << '"' << '\n';
1621	cerr << "should be: " << Signature << '\n';
1622	exit(1);
1623	}
1624
1625	if ( Data_From_STDIN )
1626	cin.read( (char *)sign, 1);
1627	else
1628	inputfile.read( (char *)sign, 1);
1629
1630	}
1631
1632	// huschel ends here
1633
1634	} // end if found end of file
1635
1636	} // end if found end of run
1637
1638	if ( Data_From_STDIN )
1639	cin.read( flag, SIZE_OF_FLAGS );
1640	else
1641	inputfile.read ( flag, SIZE_OF_FLAGS );
1642
1643	} // end if else found start of run
1644	} // end big while loop
1645
1646	//!@' @#### End of program.
1647	//@'
1648
1649	//end my version
1650
1651	#ifdef __ROOT__
1652	//++
1653	// put the Event to the root file
1654	//--
1655
1656	EvtTree.Write() ;
1657	outfile.Write() ;
1658	outfile.Close() ;
1659
1660	#endif // __ROOT__
1661
1662	// close input file
1663
1664	ntcph_in += ncph_in;
1665	ntcph_out += ncph_out;
1666	log( SIGNATURE,
1667	"%d event(s), with a total of %d C.photons in and %d C.photons out \n",
1668	ntshow, ntcph_in, ntcph_out );
1669
1670	// log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
1671	// ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
1672
1673	// close files
1674
1675	log( SIGNATURE, "Closing files\n" );
1676
1677	inputfile.close();
1678
1679	#ifdef __DETAIL_TRIGGER__
1680	// Output of Trigger statistics
1681	//
1682
1683	// Trigger.PrintStat() ;
1684	#endif // __DETAIL_TRIGGER__
1685
1686	// program finished
1687
1688	log( SIGNATURE, "Done.\n");
1689
1690	return( 0 );
1691
1692	}
1693	//!@}
1694
1695	// @T \newpage
1696
1697	//!@subsection Functions definition.
1698
1699	//!-----------------------------------------------------------
1700	// @name present
1701	//
1702	// @desc Make some presentation
1703	//
1704	// @date Sat Jun 27 05:58:56 MET DST 1998
1705	//------------------------------------------------------------
1706	// @function
1707
1708	//!@{
1709	void
1710	present(void)
1711	{
1712	cout << "##################################################\n"
1713	<< SIGNATURE << '\n' << '\n'
1714	<< "Processor of the reflector output\n"
1715	<< "J C Gonzalez, Jun 1998\n"
1716	<< "##################################################\n\n"
1717	<< flush ;
1718	}
1719	//!@}
1720
1721
1722	//!-----------------------------------------------------------
1723	// @name usage
1724	//
1725	// @desc show help
1726	//
1727	// @date Tue Dec 15 16:23:30 MET 1998
1728	//------------------------------------------------------------
1729	// @function
1730
1731	//!@{
1732	void
1733	usage(void)
1734	{
1735	present();
1736	cout << "\nusage ::\n\n"
1737	<< "\t camera "
1738	<< " [ -@ paramfile ] "
1739	<< " [ -h ] "
1740	<< "\n\n or \n\n"
1741	<< "\t camera < paramfile"
1742	<< "\n\n";
1743	exit(0);
1744	}
1745	//!@}
1746
1747
1748	//!-----------------------------------------------------------
1749	// @name log
1750	//
1751	// @desc function to send log information
1752	//
1753	// @var funct Name of the caller function
1754	// @var fmt Format to be used (message)
1755	// @var ... Other information to be shown
1756	//
1757	// @date Sat Jun 27 05:58:56 MET DST 1998
1758	//------------------------------------------------------------
1759	// @function
1760
1761	//!@{
1762	void
1763	log(const char funct, char fmt, ...)
1764	{
1765	va_list args;
1766
1767	// Display the name of the function that called error
1768	printf("[%s]: ", funct);
1769
1770	// Display the remainder of the message
1771	va_start(args, fmt);
1772	vprintf(fmt, args);
1773	va_end(args);
1774	}
1775	//!@}
1776
1777
1778	//!-----------------------------------------------------------
1779	// @name error
1780	//
1781	// @desc function to send an error message, and abort the program
1782	//
1783	// @var funct Name of the caller function
1784	// @var fmt Format to be used (message)
1785	// @var ... Other information to be shown
1786	//
1787	// @date Sat Jun 27 05:58:56 MET DST 1998
1788	//------------------------------------------------------------
1789	// @function
1790
1791	//!@{
1792	void
1793	error(const char funct, char fmt, ...)
1794	{
1795	va_list args;
1796
1797	// Display the name of the function that called error
1798	fprintf(stderr, "ERROR in %s: ", funct);
1799
1800	// Display the remainder of the message
1801	va_start(args, fmt);
1802	vfprintf(stderr, fmt, args);
1803	va_end(args);
1804
1805	perror(funct);
1806
1807	abort();
1808	}
1809	//!@}
1810
1811
1812	//!-----------------------------------------------------------
1813	// @name isA
1814	//
1815	// @desc returns TRUE(FALSE), if the flag is(is not) the given
1816	//
1817	// @var s1 String to be searched
1818	// @var flag Flag to compare with string s1
1819	// @return TRUE: both strings match; FALSE: oth.
1820	//
1821	// @date Wed Jul 8 15:25:39 MET DST 1998
1822	//------------------------------------------------------------
1823	// @function
1824
1825	//!@{
1826	int
1827	isA( char * s1, const char * flag ) {
1828	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
1829	}
1830	//!@}
1831
1832
1833	//!-----------------------------------------------------------
1834	// @name read_ct_file
1835	//
1836	// @desc read CT definition file
1837	//
1838	// @date Sat Jun 27 05:58:56 MET DST 1998
1839	//------------------------------------------------------------
1840	// @function
1841
1842	//!@{
1843	void
1844	read_ct_file(void)
1845	{
1846	char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
1847	char token[ITEM_MAX_LENGTH]; //@< a single token
1848	int i, j; //@< dummy counters
1849
1850	log( "read_ct_file", "start.\n" );
1851
1852	ifstream ctin ( ct_filename );
1853
1854	if ( ctin.bad() )
1855	error( "read_ct_file",
1856	"Cannot open CT def. file: %s\n", ct_filename );
1857
1858	// loop till the "end" directive is reached
1859
1860	while (!ctin.eof()) {
1861
1862	// get line from stdin
1863
1864	ctin.getline(line, LINE_MAX_LENGTH);
1865
1866	// look for each item at the beginning of the line
1867
1868	for (i=0; i<=define_mirrors; i++)
1869	if (strstr(line, CT_ITEM_NAMES[i]) == line)
1870	break;
1871
1872	// if it is not a valid line, just ignore it
1873
1874	if (i == define_mirrors+1)
1875	continue;
1876
1877	// case block for each directive
1878
1879	switch ( i ) {
1880
1881	case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
1882
1883	// get focal distance
1884
1885	sscanf(line, "%s %d", token, &ct_Type);
1886
1887	log( "read_ct_file", "<Type of Telescope>: %s\n",
1888	((ct_Type==0) ? "CT1" : "MAGIC") );
1889
1890	break;
1891
1892	case focal_distance: // <focal distance> [cm]
1893
1894	// get focal distance
1895
1896	sscanf(line, "%s %f", token, &ct_Focal_mean);
1897
1898	log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
1899
1900	break;
1901
1902	case focal_std: // s(focal distance) [cm]
1903
1904	// get focal distance
1905
1906	sscanf(line, "%s %f", token, &ct_Focal_std);
1907
1908	log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
1909
1910	break;
1911
1912	case point_spread: // <point spread> [cm]
1913
1914	// get point spread
1915
1916	sscanf(line, "%s %f", token, &ct_PSpread_mean);
1917
1918	log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
1919
1920	break;
1921
1922	case point_std: // s(point spread) [cm]
1923
1924	// get point spread
1925
1926	sscanf(line, "%s %f", token, &ct_PSpread_std);
1927
1928	log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
1929
1930	break;
1931
1932	case adjustment_dev: // s(adjustment_dev) [cm]
1933
1934	// get point spread
1935
1936	sscanf(line, "%s %f", token, &ct_Adjustment_std);
1937
1938	log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
1939
1940	break;
1941
1942	case black_spot: // radius of the black spot in the center [cm]
1943
1944	// get black spot radius
1945
1946	sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
1947
1948	log( "read_ct_file", "Radius of the black spots: %f cm\n",
1949	ct_BlackSpot_rad);
1950
1951	break;
1952
1953	case r_mirror: // radius of the mirrors [cm]
1954
1955	// get radius of mirror
1956
1957	sscanf(line, "%s %f", token, &ct_RMirror);
1958
1959	log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
1960
1961	break;
1962
1963	case n_mirrors: // number of mirrors
1964
1965	// get the name of the output_file from the line
1966
1967	sscanf(line, "%s %d", token, &ct_NMirrors);
1968
1969	log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
1970
1971	break;
1972
1973	case camera_width: // camera width [cm]
1974
1975	// get the name of the ct_file from the line
1976
1977	sscanf(line, "%s %f", token, &ct_CameraWidth);
1978
1979	log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
1980
1981	break;
1982
1983	case n_pixels: // number of pixels
1984
1985	// get the name of the output_file from the line
1986
1987	sscanf(line, "%s %d", token, &ct_NPixels);
1988
1989	log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
1990
1991	break;
1992
1993	case n_centralpixels: // number of central pixels
1994
1995	// get the name of the output_file from the line
1996
1997	sscanf(line, "%s %d", token, &ct_NCentralPixels);
1998
1999	log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
2000
2001	break;
2002
2003	case n_gappixels: // number of gap pixels
2004
2005	// get the name of the output_file from the line
2006
2007	sscanf(line, "%s %d", token, &ct_NGapPixels);
2008
2009	log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
2010
2011	break;
2012
2013	case pixel_width: // pixel width [cm]
2014
2015	// get the name of the ct_file from the line
2016
2017	sscanf(line, "%s %f", token, &ct_PixelWidth);
2018
2019	ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
2020	ct_PixelWidth_corner_2_corner_half =
2021	ct_PixelWidth_corner_2_corner * 0.50;
2022	ct_Apot = ct_PixelWidth / 2;
2023	ct_2Apot = ct_Apot * 2.0;
2024
2025	log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2026
2027	break;
2028
2029	case define_mirrors: // read table with the parameters of the mirrors
2030
2031	log( "read_ct_file", "Table of mirrors data:\n" );
2032
2033	// check whether the number of mirrors was already set
2034
2035	if ( ct_NMirrors == 0 )
2036	error( "read_ct_file", "NMirrors was not set.\n" );
2037
2038	// allocate memory for paths list
2039
2040	log( "read_ct_file", "Allocating memory for ct_data\n" );
2041
2042	ct_data = new float*[ct_NMirrors];
2043
2044	for (i=0; i<ct_NMirrors; i++)
2045	ct_data[i] = new float[CT_NDATA];
2046
2047	// read data
2048
2049	log( "read_ct_file", "Reading mirrors data...\n" );
2050
2051	for (i=0; i<ct_NMirrors; i++)
2052	for (j=0; j<CT_NDATA; j++)
2053	ctin >> ct_data[i][j];
2054
2055	break;
2056
2057	} // switch ( i )
2058
2059	} // end while
2060
2061	// end
2062
2063	log( "read_ct_file", "done.\n" );
2064
2065	return;
2066	}
2067	//!@}
2068
2069
2070	//!-----------------------------------------------------------
2071	// @name read_pixels
2072	//
2073	// @desc read pixels data
2074	//
2075	// @date Fri Mar 12 16:33:34 MET 1999
2076	//------------------------------------------------------------
2077	// @function
2078
2079	//!@{
2080	void
2081	read_pixels(struct camera *pcam)
2082	{
2083	ifstream qefile;
2084	char line[LINE_MAX_LENGTH];
2085	int n, i, j, k;
2086	float qe;
2087
2088	//------------------------------------------------------------
2089	// first, pixels' coordinates
2090
2091	pcam->inumpixels = ct_NPixels;
2092	pcam->inumcentralpixels = ct_NCentralPixels;
2093	pcam->inumgappixels = ct_NGapPixels;
2094	pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
2095	pcam->dpixdiameter_cm = ct_PixelWidth;
2096
2097	// initialize pixel numbers
2098
2099	for ( i=0; i<PIX_ARRAY_SIDE; ++i )
2100	for ( j=0; j<PIX_ARRAY_SIDE; ++j )
2101	pixels[i][j][PIXNUM] = -1;
2102
2103	pixary = new float* [2*ct_NCentralPixels];
2104	for ( i=0; i<2*ct_NCentralPixels; ++i )
2105	pixary[i] = new float[2];
2106
2107	pixneig = new int* [ct_NCentralPixels];
2108	for ( i=0; i<ct_NCentralPixels; ++i ) {
2109	pixneig[i] = new int[6];
2110	for ( j=0; j<6; ++j )
2111	pixneig[i][j] = -1;
2112	}
2113
2114	npixneig = new int[ct_NCentralPixels];
2115	for ( i=0; i<ct_NCentralPixels; ++i )
2116	npixneig[i] = 0;
2117
2118	// generate all coordinates
2119
2120	igen_pixel_coordinates(pcam);
2121
2122	// transfer coordinates to the working arrays for
2123	// the central pixels
2124
2125	for(k=0; k<ct_NCentralPixels; k++){
2126
2127	i = (int) pcam->di[k];
2128	j = (int) pcam->dj[k];
2129
2130	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k;
2131	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k];
2132	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k];
2133
2134	pixary[k][0] = pcam->dxc[k];
2135	pixary[k][1] = pcam->dyc[k];
2136	}
2137
2138	// calculate tables of neighbours
2139
2140	#ifdef __DEBUG__
2141	for ( n=0 ; n<ct_NPixels ; ++n ) {
2142	cout << "Para el pixel " << n << ": ";
2143	for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2144	if ( pixels_are_neig(n,i) == TRUE ) {
2145	pixneig[n][npixneig[n]] = i;
2146	pixneig[i][npixneig[i]] = n;
2147	cout << i << ' ';
2148	++npixneig[n];
2149	++npixneig[i];
2150	}
2151	}
2152	cout << endl << flush;
2153	}
2154	#else // ! __DEBUG__
2155	for ( n=0 ; n<ct_NCentralPixels ; ++n )
2156	for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2157	if ( pixels_are_neig(n,i) == TRUE ) {
2158	pixneig[n][npixneig[n]] = i;
2159	pixneig[i][npixneig[i]] = n;
2160	++npixneig[n];
2161	++npixneig[i];
2162	}
2163	#endif // ! __DEBUG__
2164
2165	#ifdef __DEBUG__
2166	for ( n=0 ; n<ct_NPixels ; ++n ) {
2167	cout << n << ':';
2168	for ( j=0; j<npixneig[n]; ++j)
2169	cout << ' ' << pixneig[n][j];
2170	cout << endl << flush;
2171	}
2172	#endif // __DEBUG__
2173
2174	//------------------------------------------------------------
2175	// second, pixels' QE
2176
2177	// try to open the file
2178
2179	log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE);
2180
2181	qefile.open( QE_FILE );
2182
2183	// if it is wrong or does not exist, exit
2184
2185	if ( qefile.bad() )
2186	error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2187
2188	// read file
2189
2190	log("read_pixels", "Reading data . . .\n");
2191
2192	i=-1;
2193
2194	while ( ! qefile.eof() ) {
2195
2196	// get line from the file
2197
2198	qefile.getline(line, LINE_MAX_LENGTH);
2199
2200	// skip if comment
2201
2202	if ( *line == '#' )
2203	continue;
2204
2205	// if it is the first valid value, it is the number of QE data points
2206
2207	if ( i < 0 ) {
2208
2209	// get the number of datapoints
2210
2211	sscanf(line, "%d", &pointsQE);
2212
2213	// allocate memory for the table of QEs
2214
2215	QE = new float ** [ct_NPixels];
2216
2217	for ( i=0; i<ct_NPixels; ++i ) {
2218	QE[i] = new float * [2];
2219	QE[i][0] = new float[pointsQE];
2220	QE[i][1] = new float[pointsQE];
2221	}
2222
2223	QElambda = new float [pointsQE];
2224
2225	for ( i=0; i<pointsQE; ++i ) {
2226	qefile.getline(line, LINE_MAX_LENGTH);
2227	sscanf(line, "%f", &QElambda[i]);
2228	}
2229
2230	i=0;
2231
2232	continue;
2233	}
2234
2235	// get the values (num-pixel, num-datapoint, QE-value)
2236
2237	sscanf(line, "%d %d %f", &i, &j, &qe);
2238
2239	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2240	((j-1) < pointsQE) && ((j-1) > -1) ) {
2241	QE[i-1][0][j-1] = QElambda[j-1];
2242	QE[i-1][1][j-1] = qe;
2243	}
2244
2245	}
2246
2247	// close file
2248
2249	qefile.close();
2250
2251	// end
2252
2253	log("read_pixels", "Done.\n");
2254
2255	}
2256	//!@}
2257
2258
2259	//!-----------------------------------------------------------
2260	// @name pixels_are_neig
2261	//
2262	// @desc check whether two pixels are neighbours
2263	//
2264	// @var pix1 Number of the first pixel
2265	// @var pix2 Number of the second pixel
2266	// @return TRUE: both pixels are neighbours; FALSE: oth.
2267	//
2268	// @date Wed Sep 9 17:58:37 MET DST 1998
2269	//------------------------------------------------------------
2270	// @function
2271
2272	//!@{
2273	int
2274	pixels_are_neig(int pix1, int pix2)
2275	{
2276	if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2277	SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2278	> ct_PixelWidth_corner_2_corner )
2279	return ( FALSE );
2280	else
2281	return ( TRUE );
2282	}
2283	//!@}
2284
2285	//!-----------------------------------------------------------
2286	// @name igen_pixel_coordinates
2287	//
2288	// @desc generate the pixel center coordinates
2289	//
2290	// @var *pcam structure camera containing all the
2291	// camera information
2292	// @return total number of pixels
2293	//
2294	// DP
2295	//
2296	// @date Thu Oct 14 10:41:03 CEST 1999
2297	//------------------------------------------------------------
2298	// @function
2299
2300	//!@{
2301	/****** igen_pixel_coordinates() *******************************/
2302
2303	int igen_pixel_coordinates(struct camera *pcam) {
2304	/* generate pixel coordinates, return value is number of pixels */
2305
2306	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2307	float fsegment_fract;
2308	double dtsize;
2309	double dhsize;
2310	double dpsize;
2311	double dxfirst_pix;
2312	double dyfirst_pix;
2313	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
2314	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
2315
2316
2317	double dstartx, dstarty; /* for the gap pixels and outer pixels */
2318	int j, nrow;
2319
2320	dpsize = pcam->dpixdiameter_cm;
2321	dtsize = dpsize * sqrt(3.) / 2.;
2322	dhsize = dpsize / 2.;
2323
2324	/* Loop over central pixels to generate co-ordinates */
2325
2326	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2327
2328	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2329
2330	pcam->dpixsizefactor[ipixno] = 1.;
2331
2332	in = 0;
2333
2334	i = 0;
2335	itot_inside_ring = 0;
2336	iring_no = 0;
2337
2338	/* Calculate the number of pixels out to and including the ring containing pixel number */
2339	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2340
2341	while (itot_inside_ring == 0){
2342
2343	iN = 3(i(i+1)) + 1;
2344
2345	if (ipixno <= iN){
2346	iring_no = i;
2347	itot_inside_ring = iN;
2348	}
2349
2350	i++;
2351	}
2352
2353
2354	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2355
2356	ipix_in_ring = 0;
2357	for (i = 0; i < iring_no; ++i){
2358
2359	ipix_in_ring = ipix_in_ring + 6;
2360	}
2361
2362	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2363	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2364	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2365
2366	isegment = 0;
2367	fsegment_fract = 0.;
2368	if (iring_no > 0) {
2369
2370	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2371
2372	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2373
2374	}
2375
2376	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2377	/* pixel width (flat to flat) dpsize. */
2378
2379	dxfirst_pix = dpsize*iring_no;
2380	dyfirst_pix = 0.;
2381
2382	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2383
2384	ddxseg1 = - dhsize*iring_no;
2385	ddyseg1 = dtsize*iring_no;
2386	ddxseg2 = -dpsize*iring_no;
2387	ddyseg2 = 0.;
2388	ddxseg3 = ddxseg1;
2389	ddyseg3 = -ddyseg1;
2390	ddxseg4 = -ddxseg1;
2391	ddyseg4 = -ddyseg1;
2392	ddxseg5 = -ddxseg2;
2393	ddyseg5 = 0.;
2394	ddxseg6 = -ddxseg1;
2395	ddyseg6 = ddyseg1;
2396
2397	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2398	/* anti-clockwise around the ring by adding the segment to segment vectors. */
2399
2400	switch (isegment) {
2401
2402	case 0:
2403
2404	pcam->dxc[ipixno-1] = 0.;
2405	pcam->dyc[ipixno-1] = 0.;
2406
2407	case 1:
2408	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2409	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2410
2411	break;
2412
2413	case 2:
2414
2415	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2416	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2417
2418	break;
2419
2420	case 3:
2421
2422	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
2423	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
2424
2425	break;
2426
2427	case 4:
2428
2429	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
2430	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
2431
2432	break;
2433
2434	case 5:
2435
2436	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
2437	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2438
2439	break;
2440
2441	case 6:
2442
2443	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2444	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2445
2446	break;
2447
2448	default:
2449
2450	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2451	return(0);
2452
2453	} /* end switch */
2454
2455	} /* end for */
2456
2457	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2458	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2459
2460	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2461
2462	j = pcam->inumcentralpixels;
2463
2464	for(i=0; i<pcam->inumgappixels; i=i+6){
2465	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2466	pcam->dyc[j + i ] = dstarty;
2467	pcam->dpixsizefactor[j + i] = 1.;
2468	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2469	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2470	pcam->dpixsizefactor[j + i + 1] = 1.;
2471	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2472	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2473	pcam->dpixsizefactor[j + i+ 2] = 1.;
2474	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2475	pcam->dyc[j + i + 3] = dstarty;
2476	pcam->dpixsizefactor[j + i+ 3] = 1.;
2477	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2478	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2479	pcam->dpixsizefactor[j + i+ 4] = 1.;
2480	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2481	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2482	pcam->dpixsizefactor[j + i + 5] = 1.;
2483	} /* end for */
2484	} /* end if */
2485
2486	/* generate positions of the outer pixels */
2487
2488	if( pcam->inumbigpixels > 0 ){
2489
2490	j = pcam->inumcentralpixels + pcam->inumgappixels;
2491
2492	for(i=0; i<pcam->inumbigpixels; i++){
2493	pcam->dpixsizefactor[j + i] = 2.;
2494	}
2495
2496	in = 0;
2497
2498	nrow = (int) ceil(dstartx / 2. / dpsize);
2499
2500	while(in < pcam->inumbigpixels){
2501
2502	pcam->dxc[j + in] = dstartx + dpsize;
2503	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2504	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2505	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2506	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2507	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2508	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2509	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2510	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2511	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2512	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2513	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2514	for(i=1; i<nrow; i++){
2515	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2516	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2517	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2518	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2519	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2520	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2521	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2522	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2523	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
2524	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
2525	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
2526	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
2527	}
2528	in = in + 6 * nrow;
2529	dstartx = dstartx + 2. * dpsize;
2530	nrow = nrow + 1;
2531	} /* end while */
2532
2533	} /* end if */
2534
2535	/* generate the ij coordinates */
2536
2537	for(i=0; i<pcam->inumpixels; i++){
2538	pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize;
2539	pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60;
2540
2541	// fprintf(stdout, "%d %f %f %f %f %f\n",
2542	// i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i],
2543	// pcam->dpixsizefactor[i]);
2544
2545	}
2546
2547	return(pcam->inumpixels);
2548
2549	}
2550	//!@}
2551
2552	//!-----------------------------------------------------------
2553	// @name bpoint_is_in_pix
2554	//
2555	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
2556	//
2557	// @var *pcam structure camera containing all the
2558	// camera information
2559	// @var dx, dy point coordinates in centimeters
2560	// @var ipixnum pixel number (starting at 0)
2561	// @return TRUE if the point is inside the pixel, FALSE otherwise
2562	//
2563	// DP
2564	//
2565	// @date Thu Oct 14 16:59:04 CEST 1999
2566	//------------------------------------------------------------
2567	// @function
2568
2569	//!@{
2570
2571	/****** bpoint_is_in_pix() *************************************/
2572
2573	int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
2574	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
2575	/* the pixel is assumed to be a "closed set" */
2576
2577	double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
2578	double c, xx, yy; /* auxiliary variable */
2579
2580	b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
2581	a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
2582	c = 1. - 1./sqrt(3.);
2583	if((ipixnum < 0)\|\|(ipixnum >= pcam->inumpixels)){
2584	fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
2585	fprintf(stderr, "Exiting.\n");
2586	exit(203);
2587	}
2588	xx = dx - pcam->dxc[ipixnum];
2589	yy = dy - pcam->dyc[ipixnum];
2590
2591	if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) \|\|
2592	((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
2593	return(TRUE); /* point is inside */
2594	}
2595	else{
2596	return(FALSE); /* point is outside */
2597	}
2598	}
2599
2600	//!@}
2601
2602	//------------------------------------------------------------
2603	// @name dist_r_P
2604	//
2605	// @desc distance straight line r - point P
2606	//
2607	// @date Sat Jun 27 05:58:56 MET DST 1998
2608	// @function @code
2609	//------------------------------------------------------------
2610	// dist_r_P
2611	//
2612	// distance straight line r - point P
2613	//------------------------------------------------------------
2614
2615	float
2616	dist_r_P(float a, float b, float c,
2617	float l, float m, float n,
2618	float x, float y, float z)
2619	{
2620	return (
2621	sqrt((SQR((a-x)m-(b-y)l) +
2622	SQR((b-y)n-(c-z)m) +
2623	SQR((c-z)l-(a-x)n))/
2624	(SQR(l)+SQR(m)+SQR(n))
2625	)
2626	);
2627	}
2628	// @endcode
2629
2630
2631	//=------------------------------------------------------------
2632	//!@subsection Log of this file.
2633
2634	//!@{
2635	//
2636	// $Log: not supported by cvs2svn $
2637	// Revision 1.1.1.1 2000/02/08 15:13:44 harald
2638	// this is just the starting point of the further development of the camera
2639	// project to simulate MonteCarloShower for MAGIC.
2640	// The TimeCam is using also the information of the arrival time of the
2641	// cerenkov photons (or the photoelectrons).
2642	// This is the status of the program as presented on the general MAGIC meeting
2643	// in Barcelona!
2644	// Each one is invited to collaborate!!!
2645	//
2646	//
2647	// Revision 1.1.1.1 1999/11/05 11:59:31 harald
2648	// This the starting point for CVS controlled further developments of the
2649	// camera program. The program was originally written by Jose Carlos.
2650	// But here you can find a "rootified" version to the program. This means
2651	// that there is no hbook stuff in it now. Also the output of the
2652	// program changed to the MagicRawDataFormat.
2653	//
2654	// The "rootification" was done by Dirk Petry and Harald Kornmayer.
2655	//
2656	// In the following you can see the README file of that version:
2657	//
2658	// ==================================================
2659	//
2660	// Fri Oct 22 1999 D.P.
2661	//
2662	// The MAGIC Monte Carlo System
2663	//
2664	// Camera Simulation Programme
2665	// ---------------------------
2666	//
2667	// 1) Description
2668	//
2669	// This version is the result of the fusion of H.K.'s
2670	// root_camera which is described below (section 2)
2671	// and another version by D.P. which had a few additional
2672	// useful features.
2673	//
2674	// The version compiles under Linux with ROOT 2.22 installed
2675	// (variable ROOTSYS has to be set).
2676	//
2677	// Compile as before simply using "make" in the root_camera
2678	// directory.
2679	//
2680	// All features of H.K.'s root_camera were retained.
2681	//
2682	// Additional features of this version are:
2683	//
2684	// a) HBOOK is no longer used and all references are removed.
2685	//
2686	// b) Instead of HBOOK, the user is given now the possibility of
2687	// having Diagnostic data in ROOT format as a complement
2688	// to the ROOT Raw data.
2689	//
2690	// This data is written to the file which is determined by
2691	// the new input parameter "diag_file" in the camera parameter
2692	// file.
2693	//
2694	// All source code file belonging to this part have filenames
2695	// starting with "MDiag".
2696	//
2697	// The user can read the output file using the following commands
2698	// in an interactive ROOT session:
2699	//
2700	// root [0] .L MDiag.so
2701	// root [1] new TFile("diag.root");
2702	// root [2] new TTreeViewer("T");
2703	//
2704	// This brings up a viewer from which all variables of the
2705	// TTree can be accessed and histogrammed. This example
2706	// assumes that you have named the file "diag.root", that
2707	// you are using ROOT version 2.22 or later and that you have
2708	// the shared object library "MDiag.so" which is produced
2709	// by the Makefile along with the executable "camera".
2710	//
2711	// ! The contents of the so-called diag file is not yet fixed.
2712	// ! At the moment it is what J.C.G. used to put into the HBOOK
2713	// ! ntuple. In future versions the moments calculation can be
2714	// ! removed and the parameter list be modified correspondingly.
2715	//
2716	// c) Now concatenated reflector files can be read. This is useful
2717	// if you have run the reflector with different parameters but
2718	// you want to continue the analysis with all reflector data
2719	// going into ONE ROOT outputfile.
2720	//
2721	// The previous camera version contained a bug which made reading
2722	// of two or more concatenated reflector files impossible.
2723	//
2724	// d) The reflector output format was changed. It is now version
2725	// 0.4 .
2726	// The change solely consists in a shortening of the flag
2727	// definition in the file
2728	//
2729	// include-MC/MCCphoton.hxx
2730	//
2731	// ! IF YOU WANT TO READ REFLECTOR FORMAT 0.3, you can easily
2732	// ! do so by recompiling camera with the previous version of
2733	// ! include-MC/MCCphoton.hxx.
2734	//
2735	// The change was necessary for saving space and better
2736	// debugging. From now on, this format can be frozen.
2737	//
2738	// ! For producing reflector output in the new format, you
2739	// ! of course have to recompile your reflector with the
2740	// ! new include-MC/MCCphoton.hxx .
2741	//
2742	// e) A first version of the pixelization with the larger
2743	// outer pixels is implemented. THIS IS NOT YET FULLY
2744	// TESTED, but first rough tests show that it works
2745	// at least to a good approximation.
2746	//
2747	// The present version implements the camera outline
2748	// with 18 "gap-pixels" and 595 pixels in total as
2749	// shown in
2750	//
2751	// http://sarastro.ifae.es/internal/home/hardware/camera/numbering.ps
2752	//
2753	// This change involved
2754	//
2755	// (i) The file pixels.dat is no longer needed. Instead
2756	// the coordinates are generated by the program itself
2757	// (takes maybe 1 second). In the file
2758	//
2759	// pixel-coords.txt
2760	//
2761	// in the same directory as this README, you find a list
2762	// of the coordinates generated by this new routine. It
2763	// has the format
2764	//
2765	// number i j x y size-factor
2766	//
2767	// where i and j are J.C.G.'s so called biaxis hexagonal
2768	// coordinates (for internal use) and x and y are the
2769	// coordinates of the pixel centers in the standard camera
2770	// coordinate system in units of centimeters. The value
2771	// of "size-factor" determines the linear size of the pixel
2772	// relative to the central pixels.
2773	//
2774	// (ii) The magic.def file has two additional parameters
2775	// which give the number of central pixels and the
2776	// number of gap pixels
2777	//
2778	// (iii) In camera.h and camera.cxx several changes were
2779	// necessary, among them the introduction of several
2780	// new functions
2781	//
2782	// The newly suggested outline with asymmetric Winston cones
2783	// will be implemented in a later version.
2784	//
2785	// f) phe files can no longer be read since this contradicts
2786	// our philosophy that the analysis should be done with other
2787	// programs like e.g. EVITA and not with "camera" itself.
2788	// This possibility was removed.
2789	//
2790	// g) ROOT is no longer invoked with an interactive interface.
2791	// In this way, camera can better be run as a batch program and
2792	// it uses less memory.
2793	//
2794	// h) small changes concerning the variable "t_chan" were necessary in
2795	// order to avoid segmentation faults: The variable is used as an
2796	// index and it went sometimes outside the limits when camera
2797	// was reading proton data. This is because the reflector files
2798	// don't contain the photons in a chronological order and also
2799	// the timespread can be considerably longer that the foreseen
2800	// digitisation timespan. Please see the source code of camera.cxx
2801	// round about line 1090.
2802	//
2803	// j) several unused variables were removed, a few warning messages
2804	// occur when you compile camera.cxx but these can be ignored at
2805	// the moment.
2806	//
2807	// In general the program is of course not finished. It still needs
2808	// debugging, proper trigger simulation, simulation of the asymmetric
2809	// version of the outer pixels, proper NSB simulation, adaption of
2810	// the diag "ntuple" contents to our need and others small improvements.
2811	//
2812	// In the directory rfl-files there is now a file in reflector format 0.4
2813	// containing a single event produced by the starfiled adder. It has
2814	// a duration of 30 ns and represents the region around the Crab Nebula.
2815	// Using the enclosed input parameter file, camera should process this
2816	// file without problems.
2817	//
2818	// 2) The README for the previous version of root_camera
2819	//
2820	// README for a preliminary version of the
2821	// root_camera program.
2822	//
2823	// root_camera is based on the program "camera"of Jose Carlos
2824	// Gonzalez. It was changed in the way that only the pixelisation
2825	// and the distibution of the phe to the FADCs works in a
2826	// first version.
2827	//
2828	// Using the #undef command most possibilities of the orignal
2829	// program are switched of.
2830	//
2831	// The new parts are signed by
2832	//
2833	// - ROOT or __ROOT__
2834	// nearly all important codelines for ROOT output are enclosed
2835	// in structures like
2836	// #ifdef __ROOT__
2837	//
2838	// code
2839	//
2840	// #endif // __ROOT__
2841	//
2842	// In same case the new lines are signed by a comment with the word
2843	// ROOT in it.
2844	//
2845	// For timing of the pulse some variable names are changed.
2846	// (t0, t1, t --> t_ini, t_fin, t_1st, t_chan,...)
2847	// Look also for this changes.
2848	//
2849	// For the new root-file is also a change in readparm-files
2850	//
2851	//
2852	// - __DETAIL_TRIGGER__
2853	//
2854	// This is for the implementation of the current work on trigger
2855	// studies. Because the class MTrigger is not well documented it
2856	// isn´t a part of this tar file. Only a dummy File exists.
2857	//
2858	//
2859	//
2860	// With all files in the archive, the root_camera program should run.
2861	//
2862	// A reflector file is in the directory rfl-files
2863	//
2864	// ==================================================
2865	//
2866	// From now on, use CVS for development!!!!
2867	//
2868	//
2869	//
2870	// Revision 1.3 1999/10/22 15:01:28 petry
2871	// version sent to H.K. and N.M. on Fri Oct 22 1999
2872	//
2873	// Revision 1.2 1999/10/22 09:44:23 petry
2874	// first synthesized version which compiles and runs without crashing;
2875	//
2876	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
2877	// first synthesised version
2878	//
2879	// Revision 1.13 1999/03/15 14:59:05 gonzalez
2880	// camera-1_1
2881	//
2882	// Revision 1.12 1999/03/02 09:56:10 gonzalez
2883	// * empty log message *
2884	//
2885	//
2886	//!@}
2887
2888	//=EOF

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

Download in other formats: