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

Last change on this file since 1674 was 1674, checked in by blanch, 22 years ago
*** empty log message ***
File size: 114.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: camera.cxx,v $
23// $Revision: 1.47 $
24// $Author: blanch $
25// $Date: 2002-12-10 17:19:31 $
26//
27////////////////////////////////////////////////////////////////////////
28// @tableofcontents @coverpage
29
30//=-----------------------------------------------------------
31//!@section Source code of |camera.cxx|.
32
33/*!@"
34
35 In this section we show the (commented) code of the program for the
36 read-out of the output files generated by the simulator of the
37 reflector, |reflector 0.3|.
38
39 @"*/
40
41//=-----------------------------------------------------------
42//!@subsection Includes and Global variables definition.
43
44//!@{
45
46// includes for ROOT
47// BEWARE: the order matters!
48
49#include "TROOT.h"
50
51#include "TRandom.h"
52#include "TApplication.h"
53
54#include "TFile.h"
55#include "TTree.h"
56#include "TBranch.h"
57#include "TCanvas.h"
58
59#include "TArrayC.h"
60
61#include "MTrigger.hxx"
62#include "MFadc.hxx"
63#include "MLons.hxx"
64
65#include "MRawRunHeader.h"
66#include "MRawEvtData.h"
67#include "MRawEvtHeader.h"
68#include "MRawCrateData.h"
69#include "MMcEvt.hxx"
70#include "MMcTrig.hxx"
71#include "MMcRunHeader.hxx"
72#include "MMcTrigHeader.hxx"
73#include "MMcFadcHeader.hxx"
74#include "MGeomCamMagic.h"
75#include "MGeomPix.h"
76
77/*!@"
78
79 All the defines are located in the file |camera.h|.
80
81 @"*/
82
83#include "camera.h"
84
85//!@}
86
87/*!@"
88
89 The following set of flags are used in time of compilation. They do
90 not affect directly the behaviour of the program at run-time
91 (though, of course, if you disconnected the option for
92 implementation of the Trigger logic, you will not be able to use any
93 trigger at all. The 'default' values mean default in the sense of
94 what you got from the server when you obtained this program.
95
96 @"*/
97
98//!@{
99
100// flag for debugging (default: OFF )
101#define __DEBUG__
102#undef __DEBUG__
103
104
105//!@}
106
107//=-----------------------------------------------------------
108//!@subsection Definition of global variables.
109
110/*!@"
111
112 Now we define some global variables with data about the telescope,
113 such as "focal distance", number of pixels/mirrors,
114 "size of the camera", and so on.
115
116 @"*/
117
118/*!@"
119
120 Depending on the telescope we are using (CT1 or MAGIC), the
121 information stored in the definition file is different.
122 The variable |ct_Type| has the value 0 when we use
123 CT1, and 1 when we use MAGIC.
124
125 @"*/
126
127//!@{
128static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC
129//!@}
130
131/*!@"
132
133 And this is the information about the whole telescope.
134
135 @"*/
136
137//!@{
138
139// parameters of the CT (from the CT definition file)
140
141////@: Focal distances [cm]
142//static float *ct_Focal;
143
144//@: Mean Focal distances [cm]
145static float ct_Focal_mean;
146
147//@: STDev. Focal distances [cm]
148static float ct_Focal_std;
149
150//@: Mean Point Spread function [cm]
151static float ct_PSpread_mean;
152
153//@: STDev. Point Spread function [cm]
154static float ct_PSpread_std;
155
156//@: STDev. Adjustmente deviation [cm]
157static float ct_Adjustment_std;
158
159//@: Radius of the Black Spot in mirror [cm]
160static float ct_BlackSpot_rad;
161
162//@: Radius of one mirror [cm]
163static float ct_RMirror;
164
165//@: Camera width [cm]
166static float ct_CameraWidth;
167
168//@: Pixel width [cm]
169static float ct_PixelWidth;
170
171//@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)
172static float ct_PixelWidth_corner_2_corner;
173
174//@: Number of mirrors
175static int ct_NMirrors = 0;
176
177//@: Number of pixels
178static int ct_NPixels;
179
180//@: Number of pixels
181static int ct_NCentralPixels;
182
183//@: Number of pixels
184static int ct_NGapPixels;
185
186//@: name of the CT definition file to use
187static char ct_filename[256];
188
189//@: list of showers to be skipped
190static int *Skip;
191
192//@: number of showers to be skipped
193static int nSkip=0;
194
195//@: flag: TRUE: data come from STDIN; FALSE: from file
196static int Data_From_STDIN = FALSE;
197
198//@: flag: TRUE: write all images to output; FALSE: only triggered showers
199static int Write_All_Images = FALSE;
200
201//@: flag: TRUE: write all data to output; FALSE: only triggered showers
202static int Write_All_Data = FALSE;
203
204static int Write_McEvt = TRUE;
205static int Write_McTrig = TRUE;
206static int Write_McFADC = TRUE;
207static int Write_RawEvt = FALSE;
208
209//@: flag: TRUE: selection on the energy
210static int Select_Energy = TRUE;
211
212//@: Lower edge of the selected energy range (in GeV)
213static float Select_Energy_le = 0.0;
214
215//@: Upper edge of the selected energy range (in GeV)
216static float Select_Energy_ue = 100000.0;
217
218//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
219static int FADC_Scan = FALSE;
220
221//@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't
222static int Trigger_Scan = FALSE;
223
224//@: flag: TRUE: loop trigger analysis over several thresholds, multiplicities and topologies; FALSE: a single trigger configuration
225static int Trigger_Loop = FALSE;
226
227//@: flag: TRUE: Different threshold for each pixel ; FALSE: same threshold for all pixels
228static int Individual_Thres_Pixel = FALSE;
229
230//@: Properties of the trigger
231static float Trigger_gate_length = 6.0;
232static float Trigger_response_ampl = 1.0;
233static float Trigger_response_fwhm = 2.0;
234static float Trigger_overlaping_time= 0.25;
235static float Trigger_noise= 0.3;
236
237//@: Properties of the FADC
238static float FADC_response_ampl = MFADC_RESPONSE_AMPLITUDE;
239static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
240static float FADC_noise = 2.0;
241
242//@: Trigger conditions for a single trigger mode
243static float qThreshold[CAMERA_PIXELS];
244static int Trigger_multiplicity = 4;
245static int Trigger_topology = 2;
246
247//@: Upper and lower edges of the trigger loop
248static float Trigger_loop_lthres = 2.0;
249static float Trigger_loop_uthres = 10.0;
250static float Trigger_loop_sthres = 1.0;
251static int Trigger_loop_lmult = 2;
252static int Trigger_loop_umult = 10;
253static int Trigger_loop_ltop = 0;
254static int Trigger_loop_utop = 2;
255
256//!@}
257
258/*!@"
259
260 The following double-pointer is a 2-dimensional table with information
261 about each pixel. The routine read_pixels will generate
262 the information for filling it using igen_pixel_coordinates().
263
264 @"*/
265
266//!@{
267// Pointer to a tables/Arrays with information about the pixels
268// and data stored on them with information about the pixels
269
270
271//@: coordinates x,y for each pixel
272static float **pixary;
273
274//@: indexes of pixels neighbours of a given one
275static int **pixneig;
276
277//@: number of neighbours a pixel have
278static int *npixneig;
279
280//@: contents of the pixels (ph.e.)
281static float *fnpix;
282
283
284//!@}
285
286/*!@"
287
288 The following double-pointer is a 2-dimensional table with the
289 Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
290 of the wavelength @$\lambda@$. The routine |read_pixels()| will read
291 also this information from the file |qe.dat|.
292
293 @"*/
294
295//!@{
296// Pointer to a table with QE, number of datapoints, and wavelengths
297
298//@: table of QE
299static float ***QE;
300
301//@: number of datapoints for the QE curve
302static int pointsQE;
303
304//@: table of QE
305static float *QElambda;
306
307//@: table of lightguide = WC;
308static float **WC;
309
310//@: number of datapoints for the WC curve
311static int pointsWC;
312
313//!@}
314
315/*!@"
316
317 The following double-pointer is a 2-dimensional table with information
318 about each mirror in the dish. The routine |read_ct_file()| will read
319 this information from the CT definition file.
320
321 @"*/
322
323//!@{
324// Pointer to a table with the following info.:
325
326static float **ct_data;
327
328/*
329 * TYPE=0 (CT1)
330 * i s rho theta x y z thetan phin xn yn zn
331 *
332 * i : number of the mirror
333 * s : arc length [cm]
334 * rho : polar rho of the position of the center of the mirror [cm]
335 * theta : polar angle of the position of the center of the mirror [cm]
336 * x : x coordinate of the center of the mirror [cm]
337 * y : y coordinate of the center of the mirror [cm]
338 * z : z coordinate of the center of the mirror [cm]
339 * thetan : polar theta angle of the direction where the mirror points to
340 * phin : polar phi angle of the direction where the mirror points to
341 * xn : xn coordinate of the normal vector in the center (normalized)
342 * yn : yn coordinate of the normal vector in the center (normalized)
343 * zn : zn coordinate of the normal vector in the center (normalized)
344 *
345 * TYPE=1 (MAGIC)
346 * i f sx sy x y z thetan phin
347 *
348 * i : number of the mirror
349 * f : focal distance of that mirror
350 * sx : curvilinear coordinate of mirror's center in X[cm]
351 * sy : curvilinear coordinate of mirror's center in X[cm]
352 * x : x coordinate of the center of the mirror [cm]
353 * y : y coordinate of the center of the mirror [cm]
354 * z : z coordinate of the center of the mirror [cm]
355 * thetan : polar theta angle of the direction where the mirror points to
356 * phin : polar phi angle of the direction where the mirror points to
357 * xn : xn coordinate of the normal vector in the center (normalized)
358 * yn : yn coordinate of the normal vector in the center (normalized)
359 * zn : zn coordinate of the normal vector in the center (normalized)
360 */
361//!@}
362
363/*!@"
364
365 We define a table into where random numbers will be stored.
366 The routines used for random number generation are provided by
367 |RANLIB| (taken from NETLIB, |www.netlib.org|), and by
368 the routine |double drand48(void)| (prototype defined in
369 |stdlib.h|) through the macro |RandomNumber| defined in
370 |camera.h|.
371
372 @"*/
373
374
375//!@}
376
377extern char FileName[];
378// switch on starfield rotation
379static int Starfield_rotate = TRUE;
380static float zenith = 0.0;//zenith angle for starfield rotation
381static float azimutal = 90.0;//azimuth angle for starfield rotation
382//=-----------------------------------------------------------
383// @subsection Main program.
384
385//!@{
386
387//++++++++++++++++++++++++++++++++++++++++
388// MAIN PROGRAM
389//----------------------------------------
390
391int main(int argc, char **argv)
392{
393
394 //!@' @#### Definition of variables.
395 //@'
396
397 char inname[256]; //@< input file name
398 char starfieldname[256]; //@< starfield input file name
399
400 char datname[256]; //@< data (ASCII) output file name
401
402 char rootname[256] ; //@< ROOT file name
403 char rootname_loop[256] ; //@< ROOT file name
404
405 char parname[256]; //@< parameters file name
406
407 char nsbpathname[256]; //@< directory with the dataabse for th NSB
408
409 char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
410 char flag_new[4]; //@< New flag for the run header in the .rfl file
411
412 int reflector_file_version; //@< vrsion of he reflector file
413
414 FILE *inputfile; //@< stream for the input file
415 ofstream datafile; //@< stream for the data file
416
417 MCEventHeader mcevth; //@< Event Header class (MC)
418 MCEventHeader_2 mcevth_2; //@< Event Header class (MC) for reflector > 0.6
419 MCCphoton cphoton; //@< Cherenkov Photon class (MC)
420
421 int inumphe; //@< number of photoelectrons in an event from showers
422 float inumphensb; //@< number of photoelectrons in an event from nsb
423
424 float arrtmin_ns; //@ arrival time of the first photoelectron
425 float arrtmax_ns; //@ arrival time of the last photoelectron
426
427 float thetaCT, phiCT; //@< parameters of a given shower
428 float thetashw, phishw; //@< parameters of a given shower
429 float coreX, coreY; //@< core position
430 float impactD; //@< impact parameter
431 float l1, m1, n1; //@< auxiliary variables
432 float l2, m2, n2; //@< auxiliary variables
433 float num, den; //@< auxiliary variables
434
435 int nshow=0; //@< partial number of shower in a given run
436 int ntshow=0; //@< total number of showers
437 int ncph=0; //@< partial number of photons in a given run
438 int ntcph=0; //@< total number of photons
439
440 int i, j, k; //@< simple counters
441
442 int addElecNoise; //@< Will we add ElecNoise?
443
444 int simulateNSB; //@< Will we simulate NSB?
445 int nphe2NSB; //@< From how many phe will we simulate NSB?
446 float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel)
447 float diffnsb_phepns[iMAXNUMPIX]; //@< diffuse NSB values for each pixel
448 //@< derived from meanNSB
449 float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS]; //@< non-diffuse nsb
450 //@< photoelectron rates
451 float nsb_phepns[iMAXNUMPIX]; //@< NSB values for each pixel
452
453 Float_t nsb_trigresp[TRIGGER_TIME_SLICES]; //@< array to write the trigger
454 //@< response from the database
455 Float_t nsb_fadcresp[(Int_t) SLICES_MFADC]; //@< array to write the fadc
456 //@< response from the database
457 Byte_t trigger_map[((Int_t)(CAMERA_PIXELS/8))+1]; //@< Pixels on when the
458 //@< camera triggers
459 Float_t fadc_elecnoise[CAMERA_PIXELS]; //@< Electronic niose for each pixel
460 Float_t fadc_pedestals[CAMERA_PIXELS]; //@< array for fadc pedestals values
461
462 float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
463 EXTWAVEBAND1,
464 EXTWAVEBAND2,
465 EXTWAVEBAND3,
466 EXTWAVEBAND4,
467 EXTWAVEBAND5
468 };
469 float zenfactor=1.0; // correction factor calculated from the extinction
470
471 int anaPixels;
472
473 float fCorrection; //@< Factor to apply to pixel values (def. 1.)
474
475 int ntrigger = 0; //@< number of triggers in the whole file
476 int btrigger = 0; //@< trigger flag
477 int ithrescount; //@< counter for loop over threshold trigger
478 float fthrescount; //@< value for loop over threshold trigger
479 int imulticount; //@< counter for loop over multiplicity trigger
480 int itopocount; //@< counter for loop over topology trigger
481 int isorttopo[3]; //@< sorting the topologies
482 int icontrigger; //@< number of trigger conditions to be analised
483
484 float fpixelthres[CAMERA_PIXELS]; //@< Threshold values
485
486 TArrayC *fadcValues; //@< the analog Fadc signal for pixels
487
488 float degTriggerZone; //@< trigger area in the camera (radius, in deg.)
489
490 int still_in_loop = FALSE;
491
492 //@< Variables to fill the McRunHeader
493 Int_t sfRaH = 5;
494 Int_t sfRaM = 34;
495 Int_t sfRaS = 32;
496 Int_t sfDeD = 22;
497 Int_t sfDeM = 00;
498 Int_t sfDeS = 52;
499 Float_t shthetamax = 0.0;
500 Float_t shthetamin = 0.0;
501 Float_t shphimax = 0.0;
502 Float_t shphimin = 0.0;
503 Float_t telestheta = 0.0;
504 Float_t telesphi = 0.0;
505 Float_t sofftheta = 0.0;
506 Float_t soffphi = 0.0;
507 UInt_t corsika = 5200 ;
508
509 MGeomCamMagic camgeom; // structure holding the camera definition
510
511 ct_NPixels=camgeom.GetNumPixels();
512 read_QE();
513 read_WC();
514
515 //!@' @#### Definition of variables for |getopt()|.
516 //@'
517
518 int ch, errflg = 0; //@< used by getopt
519
520 //!@' @#### Initialisation of some variables
521 //@'
522
523 for(i=0;i<iMAXNUMPIX;i++){
524 nsb_phepns[i]=0;
525 for(j=0;j<iNUMWAVEBANDS;j++)
526 nsbrate_phepns[i][j]=0.0; //@< Starlight rates initialised at 0.0
527 }
528 /*!@'
529
530 @#### Beginning of the program.
531
532 We start with the main program. First we (could) make some
533 presentation, and follows the reading of the parameters file (now
534 from the |stdin|), the reading of the CT parameters file, and the
535 creation of the output file, where the processed data will be
536 stored.
537
538 */
539
540 //++
541 // START
542 //--
543
544 // make unbuffered output
545
546 cout.setf ( ios::stdio );
547
548 // parse command line options (see reflector.h)
549
550 parname[0] = '\0';
551
552 optarg = NULL;
553 while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
554 switch (ch) {
555 case 'f':
556 strcpy(parname, optarg);
557 break;
558 case 'h':
559 usage();
560 break;
561 default :
562 errflg++;
563 }
564
565 // show help if error
566
567 if ( errflg>0 )
568 usage();
569
570 // make some sort of presentation
571
572 present();
573
574 // read parameters file
575
576 if ( strlen(parname) < 1 )
577 readparam(NULL);
578 else
579 readparam(parname);
580
581 // read data from file or from STDIN?
582
583 Data_From_STDIN = get_data_from_stdin();
584
585 // write all images, even those without trigger?
586
587 Write_All_Images = get_write_all_events();
588
589 // write all data (i.e., ph.e.s in pixels)
590
591 Write_All_Data = get_write_all_data();
592
593 Write_McEvt = get_write_McEvt() ;
594 Write_McTrig = get_write_McTrig() ;
595 Write_McFADC = get_write_McFadc() ;
596 Write_RawEvt = get_write_RawEvt() ;
597
598 FADC_Scan = get_FADC_Scan();
599 Trigger_Scan = get_Trigger_Scan();
600
601 Individual_Thres_Pixel = get_indi_thres_pixel();
602
603 get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm);
604
605 get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);
606
607 Trigger_Loop = get_Trigger_Loop(&Trigger_loop_lthres, &Trigger_loop_uthres, &Trigger_loop_sthres, &Trigger_loop_lmult, &Trigger_loop_umult, &Trigger_loop_ltop, &Trigger_loop_utop);
608
609 icontrigger =((int)((Trigger_loop_uthres-Trigger_loop_lthres)
610 /Trigger_loop_sthres)+1)*
611 (Trigger_loop_umult-Trigger_loop_lmult+1)*
612 (Trigger_loop_utop-Trigger_loop_ltop+1);
613
614 if (!Trigger_Loop){
615 get_Trigger_Single (qThreshold, &Trigger_multiplicity, &Trigger_topology);
616 icontrigger=1;
617 }
618 else
619 get_threshold(qThreshold);
620
621 // get filenames
622
623 strcpy( inname, get_input_filename() );
624 strcpy( starfieldname, get_starfield_filename() );
625 strcpy( datname, get_data_filename() );
626 strcpy( rootname, get_root_filename() );
627 strcpy( rootname_loop, get_loop_filename() );
628 strcpy( ct_filename, get_ct_filename() );
629 strcpy( nsbpathname, get_nsb_directory() );
630
631 // get different parameters of the simulation
632
633 addElecNoise = add_elec_noise(&FADC_noise, &Trigger_noise);
634 simulateNSB = get_nsb( &meanNSB, &nphe2NSB );
635
636 // get selections on the parameters
637
638 Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
639
640 //parameters for starfield rotation
641 get_teles_axis(&zenith, &azimutal);
642 Starfield_rotate = get_starfield_rotate();
643 log(SIGNATURE,
644 "%s:\n\t%20s: %s\n\t%20s: %f\n\t%20s: %f \n",
645 "Starfield Rotate",
646 "Rotate Starfield", ONoff(Starfield_rotate),
647 "Zenith Angle",zenith,
648 "Azimutal Angle",azimutal
649 );
650
651 // log filenames information
652
653 log(SIGNATURE,
654 "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n",
655 "Filenames",
656 "In", inname,
657 "Stars", starfieldname,
658 "NSB database", nsbpathname,
659 "CT", ct_filename,
660 "Data", datname,
661 "ROOT", rootname
662 );
663
664 // log Trigger information
665
666 if (Trigger_Loop) {
667 log(SIGNATURE,
668 "%s:\n\t%20s: from %5.2f to %5.2f by %5.2f step\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s\n",
669 "Trigger Loop mode",
670 "Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres,
671 "Multiplicity",Trigger_loop_lmult,Trigger_loop_umult,
672 "Topology",Trigger_loop_ltop,Trigger_loop_utop,
673 rootname_loop);
674 }
675 else if (Individual_Thres_Pixel == FALSE){
676 log(SIGNATURE,
677 "%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n",
678 "Single Trigger mode",
679 "Threshold",qThreshold[0],
680 "Multiplicity",Trigger_multiplicity,
681 "Topology",Trigger_topology);
682 }
683 else{
684 log(SIGNATURE,
685 "%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n",
686 "Single Trigger mode",
687 "Threshold","Different for each pixel",
688 "Multiplicity",Trigger_multiplicity,
689 "Topology",Trigger_topology);
690 }
691 // log flags information
692
693 log(SIGNATURE,
694 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
695 "Flags",
696 "Data_From_STDIN", ONoff(Data_From_STDIN),
697 "Write_All_Events", ONoff(Write_All_Images),
698 "Write_All_Data", ONoff(Write_All_Data));
699
700 // log flags information
701
702 log(SIGNATURE,
703 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
704 "Root output",
705 "Write_McEvt", ONoff(Write_McEvt),
706 "Write_McTrig", ONoff(Write_McTrig),
707 "Write_McFADC", ONoff(Write_McFADC),
708 "Write_RawEvt", ONoff(Write_RawEvt));
709
710 // log parameters information
711
712 log(SIGNATURE,
713 "%s:\n\t%20s: %f\n",
714 "Parameters",
715 "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB));
716
717 // log selections
718
719 log(SIGNATURE,
720 "%s:\n\t%20s: %s (%f:%f)\n",
721 "Selections:",
722 "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
723
724 // Definition and initialization of array to save trigger statistics
725
726 int ***ntriggerloop;
727
728 ntriggerloop= new int ** [(int)((Trigger_loop_uthres-Trigger_loop_lthres)
729 /Trigger_loop_sthres)];
730 for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;fthrescount+=Trigger_loop_sthres, ithrescount++){
731 ntriggerloop[ithrescount]= new int * [Trigger_loop_umult-Trigger_loop_lmult+1];
732 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
733 ntriggerloop[ithrescount][imulticount]= new int [Trigger_loop_utop-Trigger_loop_ltop+1];
734 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
735 ntriggerloop[ithrescount][imulticount][itopocount]=0;
736 }
737 }
738 }
739
740 // We should be careful that topologies are sort from
741 // the less to the more restrictive one.
742
743 if (Trigger_loop_utop==Trigger_loop_ltop)
744 for(int is=0; is<3;is++)
745 isorttopo[is]=is;
746 else {
747 isorttopo[0]=1;
748 isorttopo[1]=0;
749 isorttopo[2]=2;
750 }
751
752 // set all random numbers seeds
753
754 setall( get_seeds(0), get_seeds(1) );
755
756 // get list of showers to evt. skip
757
758 nSkip = get_nskip_showers();
759
760 if (nSkip > 0) {
761 Skip = new int[ nSkip ];
762 get_skip_showers( Skip );
763
764 log(SIGNATURE, "There are some showers to skip:\n");
765 for (i=0; i<nSkip; ++i)
766 log(SIGNATURE, "\tshower # %d\n", Skip[i]);
767 }
768
769
770 // read parameters from the ct.def file
771
772 read_ct_file();
773
774 Int_t Lev0, Lev1;
775
776 fadcValues = new TArrayC(FADC_SLICES);
777
778 // number of pixels for parameters
779
780 anaPixels = get_ana_pixels();
781 anaPixels = (anaPixels == -1) ? camgeom.GetNumPixels() : anaPixels;
782
783 // Switched off writing TObject
784
785 MArray::Class()->IgnoreTObjectStreamer();
786 MParContainer::Class()->IgnoreTObjectStreamer();
787
788 // initialise ROOT
789
790 TROOT simple("simple", "MAGIC Telescope Monte Carlo");
791
792 // initialise instance of Trigger and FADC classes
793
794 MTrigger Trigger(Trigger_gate_length,
795 Trigger_overlaping_time,
796 Trigger_response_ampl,
797 Trigger_response_fwhm); //@< A instance of the Class MTrigger
798
799 // Set Right Discriminator threshold, taking into account trigger pixels
800
801 Trigger.CheckThreshold(&qThreshold[0]);
802
803
804 // Initialise McTrig information class if we want to save trigger informtion
805
806 MMcTrig **McTrig = NULL;
807 MMcTrigHeader **HeaderTrig = NULL;
808 MMcFadcHeader **HeaderFadc = NULL;
809
810
811 if (Write_McTrig){
812
813 McTrig = new MMcTrig * [icontrigger];
814
815 for (i=0;i<icontrigger;i++) {
816 McTrig[i] = new MMcTrig();
817 }
818
819 HeaderTrig = new MMcTrigHeader * [icontrigger];
820
821 for (i=0;i<icontrigger;i++) {
822 HeaderTrig[i] = new MMcTrigHeader();
823 }
824 }
825
826 if (Write_McFADC){
827
828 HeaderFadc = new MMcFadcHeader * [icontrigger];
829
830 for (i=0;i<icontrigger;i++) {
831 HeaderFadc[i] = new MMcFadcHeader();
832 }
833 }
834
835 MFadc fadc(FADC_response_ampl,FADC_response_fwhm) ; //@< A instance of the Class MFadc
836
837 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
838 //
839 // Set the FADC pedestals for that run
840 // Some modifications
841 // mut be done to simulate a more realistic distribution of the pedestals.
842 // This simualtion is done int the SetPedestals methode inside the
843 // class MFadc
844 // Currentlly a input_pedestal array is declared with the pedestals.
845 // Thy can also be set randomly following a flat distribution.
846 //
847 /////////////////////////////////////////////////////////////////////
848
849 Float_t input_pedestals[CAMERA_PIXELS];
850
851 for(i=0;i<CAMERA_PIXELS;i++)
852 input_pedestals[i]=10.5;
853
854 fadc.SetPedestals(input_pedestals);
855
856 // Prepare the raw data output
857
858 // Header Tree
859
860 MRawRunHeader *RunHeader= new MRawRunHeader();
861 MMcRunHeader *McRunHeader = new MMcRunHeader();
862
863 // Header branch
864
865 MRawEvtHeader **EvtHeader = NULL;
866
867 if (Write_RawEvt) {
868 EvtHeader = new MRawEvtHeader * [icontrigger];
869
870 for (i=0;i<icontrigger;i++) {
871 EvtHeader[i] = new MRawEvtHeader();
872 }
873 }
874
875 // Data branch
876
877 MRawEvtData **EvtData = NULL; // Data branch
878
879 if (Write_RawEvt) {
880 EvtData = new MRawEvtData * [icontrigger];
881
882 for (i=0;i<icontrigger;i++) {
883 EvtData[i] = new MRawEvtData();
884 EvtData[i]->Init(RunHeader); // We need thr RunHeader to read
885 // number of pixels
886 }
887 }
888
889 MMcEvt *McEvt = new MMcEvt ();
890
891 //
892 // initalize a temporal ROOT file
893 //
894
895 TFile outfile_temp ( rootname , "RECREATE" );
896
897 // create a Tree for the Header Event
898 TTree HeaderTree("RunHeaders","Headers of Run");
899
900 // define branches of Header Tree
901
902 char help[4];
903
904 HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
905 &RunHeader);
906
907 HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
908 &McRunHeader);
909
910 if(!Trigger_Loop && Write_McTrig){
911
912 HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader",
913 &HeaderTrig[0]);
914 }
915 if (Trigger_Loop && Write_McTrig){
916 for(char branchname[20],i=0;i<icontrigger;i++){
917
918 sprintf(help,"%i",i+1);
919 strcpy (branchname, "MMcTrigHeader;");
920 strcat (branchname, & help[0]);
921 strcat (branchname, ".");
922 HeaderTree.Branch(branchname,"MMcTrigHeader",
923 &HeaderTrig[i]);
924 }
925 }
926
927 if(!Trigger_Loop && Write_McFADC){
928
929 HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader",
930 &HeaderFadc[0]);
931 }
932 if (Trigger_Loop && Write_McFADC){
933 for(char branchname[20],i=0;i<icontrigger;i++){
934
935 sprintf(help,"%i",i+1);
936 strcpy (branchname, "MMcFadcHeader;");
937 strcat (branchname, & help[0]);
938 strcat (branchname, ".");
939 HeaderTree.Branch(branchname,"MMcFadcHeader",
940 &HeaderFadc[i]);
941 }
942 }
943
944 // Fill branches for MRawRunHeader
945
946 RunHeader->SetMagicNumber(kMagicNumber);
947 RunHeader->SetFormatVersion(3);
948 RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
949 RunHeader->SetRunType(256);
950 RunHeader->SetRunNumber(0);
951 RunHeader->SetNumSamples(0,FADC_SLICES);
952 RunHeader->SetNumCrates(1);
953 RunHeader->SetNumPixInCrate(CAMERA_PIXELS);
954
955
956 // Fill branches for MMcTrigHeader
957
958 if(!Trigger_Loop && Write_McTrig){
959
960 HeaderTrig[0]->SetTopology((Short_t) Trigger_topology);
961 HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity);
962 HeaderTrig[0]->SetThreshold(qThreshold);
963 HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
964 HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
965 HeaderTrig[0]->SetOverlap(Trigger_overlaping_time);
966 HeaderTrig[0]->SetGate(Trigger_gate_length);
967 HeaderTrig[0]->SetElecNoise(Trigger_noise);
968 }
969 if(Trigger_Loop && Write_McTrig){
970
971 int iconcount;
972 for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++,fthrescount+=Trigger_loop_sthres){
973 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
974 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
975 HeaderTrig[iconcount]->SetTopology((Short_t) isorttopo[itopocount+Trigger_loop_ltop]);
976 HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
977 for(int i=0;i<CAMERA_PIXELS;i++){
978 fpixelthres[i]=
979 ((Float_t)(fthrescount)>=qThreshold[i])?
980 (Float_t)(fthrescount):qThreshold[i];
981 }
982 HeaderTrig[iconcount]->SetThreshold( fpixelthres);
983 HeaderTrig[iconcount]->SetAmplitud(Trigger_response_ampl);
984 HeaderTrig[iconcount]->SetFwhm(Trigger_response_fwhm);
985 HeaderTrig[iconcount]->SetOverlap(Trigger_overlaping_time);
986 HeaderTrig[iconcount]->SetGate(Trigger_gate_length);
987 HeaderTrig[iconcount]->SetElecNoise(Trigger_noise);
988 iconcount++;
989 }
990 }
991 }
992 }
993
994 // Fill branches for MMcFadcHeader
995
996 for(i=0;i<CAMERA_PIXELS;i++){
997 fadc_elecnoise[i]=FADC_noise;
998 }
999
1000 fadc.GetPedestals(&fadc_pedestals[0]);
1001
1002 if(!Trigger_Loop && Write_McFADC){
1003
1004 HeaderFadc[0]->SetShape(0.0);
1005 HeaderFadc[0]->SetAmplitud(FADC_response_ampl);
1006 HeaderFadc[0]->SetFwhm(FADC_response_fwhm);
1007 HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],anaPixels);
1008 HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
1009 }
1010 if(Trigger_Loop && Write_McFADC){
1011 int iconcount;
1012 for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1013 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
1014 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
1015 HeaderFadc[iconcount]->SetShape(0.0);
1016 HeaderFadc[iconcount]->SetAmplitud(Trigger_response_ampl);
1017 HeaderFadc[iconcount]->SetFwhm(Trigger_response_fwhm);
1018 HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0],anaPixels);
1019 HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
1020 iconcount++;
1021 }
1022 }
1023 }
1024 }
1025
1026 // Fill the Header Tree with the current leaves of each branch
1027 // HeaderTree.Fill() ;
1028
1029
1030 // create a Tree for the Event data stream
1031 TTree EvtTree("Events","Normal Triggered Events");
1032
1033 if (Write_McEvt){
1034
1035 EvtTree.Branch("MMcEvt","MMcEvt",
1036 &McEvt);
1037 }
1038
1039 if(!Trigger_Loop){
1040
1041 if (Write_RawEvt){
1042 EvtTree.Branch("MRawEvtHeader","MRawEvtHeader",
1043 &EvtHeader[0]);
1044 EvtTree.Branch("MRawEvtData","MRawEvtData",
1045 &EvtData[0]);
1046 }
1047 if (Write_McTrig){
1048 EvtTree.Branch("MMcTrig","MMcTrig",
1049 &McTrig[0]);
1050 }
1051 }
1052 else{
1053 if (Write_McTrig){
1054 for(char branchname[10],i=0;i<icontrigger;i++){
1055
1056 sprintf(help,"%i",i+1);
1057 strcpy (branchname, "MMcTrig;");
1058 strcat (branchname, & help[0]);
1059 strcat (branchname, ".");
1060 EvtTree.Branch(branchname,"MMcTrig",
1061 &McTrig[i]);
1062 }
1063 }
1064 }
1065
1066 if (Trigger_Loop && Write_RawEvt){
1067 for(char branchname[15],i=0;i<icontrigger;i++){
1068
1069 sprintf(help,"%i",i+1);
1070 strcpy (branchname, "MRawEvtHeader;");
1071 strcat (branchname, & help[0]);
1072 strcat (branchname, ".");
1073 EvtTree.Branch(branchname,"MRawEvtHeader",
1074 &EvtHeader[i]);
1075 }
1076 for(char branchname[15],i=0;i<icontrigger;i++){
1077
1078 sprintf(help,"%i",i+1);
1079 strcpy (branchname, "MRawEvtData;");
1080 strcat (branchname, & help[0]);
1081 strcat (branchname, ".");
1082 EvtTree.Branch(branchname,"MRawEvtData",
1083 &EvtData[i]);
1084 }
1085 }
1086
1087 TApplication theAppTrigger("App", &argc, argv);
1088
1089 if(FADC_Scan){
1090 if (gROOT->IsBatch()) {
1091 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1092 // return 1;
1093 }
1094 }
1095 if(FADC_Scan){
1096 //TApplication theAppFadc("App", &argc, argv);
1097
1098 if (gROOT->IsBatch()) {
1099 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1100 // return 1;
1101 }
1102 }
1103
1104 // set trigger area (deg)
1105
1106 if ( ! get_trigger_radius( &degTriggerZone ) )
1107 degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
1108
1109 if ( ! get_correction( &fCorrection ) )
1110 fCorrection = 1.0;
1111
1112 // prepare the NSB simulation
1113
1114 // Instance of the Mlons class
1115 MLons lons(Trigger_response_ampl, Trigger_response_fwhm,
1116 MFADC_RESPONSE_AMPLITUDE, MFADC_RESPONSE_FWHM);
1117
1118 lons.SetPath(nsbpathname);
1119
1120 if( simulateNSB ){
1121
1122 //
1123 // Calculate the non-diffuse NSB photoelectron rates
1124 //
1125
1126 k = produce_nsbrates( starfieldname,
1127 &camgeom,
1128 nsbrate_phepns);
1129
1130 if (k != 0){
1131 cout << "Error when reading starfield... \nExiting.\n";
1132 exit(1);
1133 }
1134
1135 // calculate diffuse rate correcting for the pixel size
1136
1137 const double size_inner = camgeom[0].GetR();
1138
1139 for(UInt_t ui=0; ui<camgeom.GetNumPixels(); ui++){
1140 const double actual_size = camgeom[ui].GetR();
1141 diffnsb_phepns[ui] = meanNSB *
1142 actual_size*actual_size/(size_inner*size_inner);
1143 }
1144
1145 // calculate nsb rate including diffuse and starlight
1146 // we also include the extinction effect
1147 for(j=0;j<iNUMWAVEBANDS;j++){
1148 // calculate the effect of the atmospheric extinction
1149
1150 zenfactor = pow(10., -0.4 * ext[j] );
1151
1152 for(UInt_t ui=0; ui<camgeom.GetNumPixels();ui++){
1153 nsb_phepns[ui]+=diffnsb_phepns[ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
1154 }
1155 }
1156 }
1157
1158 //
1159 // Read the reflector file with the Cherenkov data
1160 //
1161
1162 // select input file
1163
1164 if ( Data_From_STDIN ) {
1165
1166 inputfile = stdin;
1167
1168 }
1169 else{
1170
1171 log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
1172 inputfile = fopen( inname, "r" );
1173 if ( inputfile == NULL )
1174 error( SIGNATURE, "Cannot open input file: %s\n", inname );
1175
1176 }
1177
1178 // get signature, and check it
1179
1180 if((reflector_file_version=check_reflector_file( inputfile ))==FALSE){
1181 exit(1);
1182 }
1183
1184 // open data file
1185
1186 log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
1187 datafile.open( datname );
1188
1189 if ( datafile.bad() )
1190 error( SIGNATURE, "Cannot open data file: %s\n", datname );
1191
1192 // initializes flag
1193
1194 strcpy( flag, " \0" );
1195
1196 // allocate space for PMTs numbers of pixels
1197
1198 fnpix = new float [ camgeom.GetNumPixels() ];
1199
1200
1201 //!@' @#### Main loop.
1202 //@'
1203
1204 // get flag
1205
1206 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1207
1208 // loop over the file
1209
1210 still_in_loop = TRUE;
1211
1212 while (
1213 ((! Data_From_STDIN) && ( !feof(inputfile) ))
1214 ||
1215 (Data_From_STDIN && still_in_loop)
1216 ) {
1217
1218 // reading .rfl files
1219 if(!isA( flag, FLAG_START_OF_RUN )){
1220
1221 // We break the main loop
1222 cout<<"Warning: Expected start of run flag, but found:"<<flag<<endl;
1223 cout<<" We break the main loop"<<endl;
1224 break;
1225 }
1226 else { // found start of run
1227
1228 fread( flag_new, 4, 1, inputfile );
1229
1230 if(!isA( flag_new, FLAG_START_OF_HEADER)){
1231
1232 // We break the main loop
1233 cout<<"Warning: Expected start of run header flag, but found:"<<flag<<endl;
1234 cout<<" We break the main loop"<<endl;
1235 break;
1236 }
1237
1238 Float_t flag_temp[SIZE_OF_MCRUNHEADER];
1239
1240 fread( flag_temp, (SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile );
1241
1242 nshow=0;
1243
1244 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1245
1246 while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
1247 fread( flag_new, 4, 1, inputfile );
1248
1249 if(!isA( flag_new, FLAG_EVENT_HEADER)){
1250
1251 // We break while events loop
1252 cout<<"Warning: Expected start of event header flag, but found:"<<flag<<endl;
1253 cout<<" We break the while events loop"<<endl;
1254 break;
1255 }
1256
1257 //
1258 // Clear Trigger and Fadc
1259 //
1260 Trigger.Reset() ;
1261 Trigger.ClearFirst();
1262 Trigger.ClearZero();
1263 fadc.Reset() ;
1264
1265 ++nshow;
1266 log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
1267
1268 // get MCEventHeader
1269 if (reflector_file_version<6)
1270 fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
1271 else
1272 fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile );
1273
1274
1275 // calculate impact parameter (shortest distance betwee the original
1276 // trajectory of the primary (assumed shower-axis) and the
1277 // direction where the telescope points to
1278 //
1279 // we use the following equation, given that the shower core position
1280 // is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
1281 // and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
1282 // and (l2,m2,n2)
1283 //
1284 // | |
1285 // | x1-x2 y1-y2 z1-z2 |
1286 // | |
1287 // + | l1 m1 n1 |
1288 // - | |
1289 // | l2 m2 n2 |
1290 // | |
1291 // dist = ------------------------------------ ( > 0 )
1292 // [ |l1 m1|2 |m1 n1|2 |n1 l1|2 ]1/2
1293 // [ | | + | | + | | ]
1294 // [ |l2 m2| |m2 n2| |n2 l2| ]
1295 //
1296 // playing a little bit, we get this reduced for in our case:
1297 //
1298 //
1299 // dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
1300 // [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
1301 // 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
1302
1303 // read the direction of the incoming shower
1304
1305 if (reflector_file_version<6){
1306 thetashw = mcevth.get_theta();
1307 phishw = mcevth.get_phi();
1308 }
1309 else{
1310 thetashw = mcevth_2.get_theta();
1311 phishw = mcevth_2.get_phi();
1312 }
1313
1314 // calculate vector for shower
1315
1316 l1 = sin(thetashw)*cos(phishw);
1317 m1 = sin(thetashw)*sin(phishw);
1318 n1 = cos(thetashw);
1319
1320 // read the deviation of the telescope with respect to the shower
1321
1322 if (reflector_file_version<6)
1323 mcevth.get_deviations ( &thetaCT, &phiCT );
1324 else
1325 mcevth_2.get_deviations ( &thetaCT, &phiCT );
1326
1327 if ( (thetaCT == 0.) && (phiCT == 0.) ) {
1328
1329 // CT was looking to the source (both lines are parallel)
1330 // therefore, we calculate the impact parameter as the distance
1331 // between the CT axis and the core position
1332
1333 if (reflector_file_version<6)
1334 mcevth.get_core(&coreX, &coreY);
1335 else
1336 mcevth_2.get_core(&coreX, &coreY);
1337
1338 impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1339
1340 } else {
1341
1342 // the shower comes off-axis
1343
1344 // obtain with this the final direction of the CT
1345
1346 thetaCT += thetashw;
1347 phiCT += phishw;
1348
1349 // calculate vector for telescope
1350
1351 l2 = sin(thetaCT)*cos(phiCT);
1352 m2 = sin(thetaCT)*sin(phiCT);
1353 n2 = cos(thetaCT);
1354
1355 num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
1356 den = (SQR(l1*m2 - l2*m1) +
1357 SQR(m1*n2 - m2*n1) +
1358 SQR(n1*l2 - n2*l1));
1359 den = sqrt(den);
1360
1361 impactD = fabs(num)/den;
1362
1363 }
1364
1365 // energy cut
1366
1367 if ( Select_Energy ) {
1368 if (reflector_file_version<6)
1369 if (( mcevth.get_energy() < Select_Energy_le ) ||
1370 ( mcevth.get_energy() > Select_Energy_ue )) {
1371 log(SIGNATURE, "select_energy: shower rejected.\n");
1372 continue;
1373 }
1374 else
1375 if (( mcevth_2.get_energy() < Select_Energy_le ) ||
1376 ( mcevth_2.get_energy() > Select_Energy_ue )) {
1377 log(SIGNATURE, "select_energy: shower rejected.\n");
1378 continue;
1379 }
1380 }
1381
1382 // Read first and last time and put inumphe to 0
1383
1384 if (reflector_file_version<6)
1385 mcevth.get_times(&arrtmin_ns,&arrtmax_ns);
1386 else
1387 mcevth_2.get_times(&arrtmin_ns,&arrtmax_ns);
1388
1389 inumphe=0;
1390
1391 // read the photons and produce the photoelectrons
1392
1393 k = produce_phes( inputfile,
1394 &camgeom,
1395 WAVEBANDBOUND1,
1396 WAVEBANDBOUND6,
1397 &Trigger, // will be changed by the function!
1398 &fadc, // will be changed by the function!
1399 &inumphe, // important for later: the size of photoe[]
1400 fnpix, // will be changed by the function!
1401 &ncph, // will be changed by the function!
1402 &arrtmin_ns, // will be changed by the function!
1403 &arrtmax_ns, // will be changed by the function!
1404 -1 ); // photons from the shower
1405
1406 if( k != 0 ){ // non-zero returnvalue means error
1407 cout << "Exiting.\n";
1408 exit(1);
1409 }
1410
1411 // NSB simulation
1412
1413 if(simulateNSB && nphe2NSB<=inumphe){
1414
1415 // Fill trigger and fadc response in the trigger class from the database
1416 for(UInt_t ui=0;ui<camgeom.GetNumPixels();ui++)
1417 if(nsb_phepns[ui]>0.0){
1418 k=lons.GetResponse(nsb_phepns[ui],0.1,
1419 & nsb_trigresp[0],& nsb_fadcresp[0]);
1420 if(k==0){
1421 cout << "Exiting.\n";
1422 exit(1);
1423 }
1424 Trigger.AddNSB(ui,nsb_trigresp);
1425 fadc.AddSignal(ui,nsb_fadcresp);
1426 }
1427 }// end if(simulateNSB && nphe2NSB<=inumphe) ...
1428
1429 log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1430 ncph, ntcph);
1431
1432 ntcph += ncph;
1433
1434 // skip it ?
1435
1436 for ( i=0; i<nSkip; ++i ) {
1437 if (Skip[i] == (nshow+ntshow)) {
1438 i = -1;
1439 break;
1440 }
1441 }
1442
1443 // if after the previous loop, the exit value of i is -1
1444 // then the shower number is in the list of showers to be
1445 // skipped
1446
1447 if (i == -1) {
1448 log(SIGNATURE, "\t\tskipped!\n");
1449 continue;
1450 }
1451
1452 cout << "Total number of phes: " << inumphe<<" + ";
1453 inumphensb=0;
1454 for (UInt_t ui=0;ui<camgeom.GetNumPixels();ui++){
1455 inumphensb+=nsb_phepns[ui]*TOTAL_TRIGGER_TIME;
1456 }
1457 cout<<inumphensb<<endl;
1458
1459 //++++++++++++++++++++++++++++++++++++++++++++++++++
1460 // at this point we have a camera full of
1461 // ph.e.s
1462 // we should first apply the trigger condition,
1463 // and if there's trigger, then clean the image,
1464 // calculate the islands statistics and the
1465 // other parameters of the image (Hillas' parameters
1466 // and so on).
1467 //--------------------------------------------------
1468
1469 // TRIGGER HERE
1470
1471 //
1472 // now the noise of the electronic
1473 // (preamps, optical transmission,..) is introduced.
1474 // This is done inside the class MTrigger by the method ElecNoise.
1475 //
1476 if (addElecNoise){
1477 Trigger.ElecNoise(Trigger_noise) ;
1478
1479 fadc.ElecNoise(FADC_noise) ;
1480 }
1481
1482 // We should simulate the AC coupling behaviour:
1483 // For the FADC it is only done for the NSB while producinr
1484 // the StarResponse database.
1485 // For the trigger is done in the Trigger.Diskriminate(), which
1486 // is called later (it should be seperated to speed up the program
1487
1488 // now a shift in the fadc signal due to the pedestlas is
1489 // intorduced
1490 // This is done inside the class MFadc by the method Pedestals
1491 fadc.Pedestals();
1492
1493 // We study several trigger conditons
1494 if(Trigger_Loop){
1495
1496 // Set to zero the flag to know if some conditon has triggered
1497 btrigger=0;
1498 // Loop over trigger threshold
1499 int iconcount;
1500 for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1501 for (i=0;i<CAMERA_PIXELS;i++)
1502 fpixelthres[i]=
1503 ((Float_t)(fthrescount)>=qThreshold[i])?
1504 (Float_t)(fthrescount):qThreshold[i];
1505 Trigger.SetThreshold(fpixelthres);
1506
1507 Trigger.Diskriminate();
1508 //
1509 // look if in all the signals in the trigger signal branch
1510 // is a possible Trigger. Therefore we habe to diskriminate all
1511 // the simulated analog signals (Method Diskriminate in class
1512 // MTrigger). We look simultanously for the moments at which
1513 // there are more than TRIGGER_MULTI pixels above the
1514 // CHANNEL_THRESHOLD.
1515 //
1516
1517 // Set trigger flags to zero
1518 Lev1=0;
1519
1520 // loop over multiplicity of trigger configuration
1521 for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1522 Trigger.SetMultiplicity(imulticount);
1523 Trigger.ClearZero();
1524
1525 Lev0=(Short_t) Trigger.ZeroLevel();
1526 if (Lev0>0 || Write_All_Images || btrigger){
1527
1528 // loop over topologies
1529 for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1530 Lev1=0;
1531
1532 if(itopocount==0 && imulticount>7) continue;
1533 //COBB if(itopocount==2 && imulticount<3) continue;
1534 // It only makes to look for a different topology
1535 // if there are 3 or more N pixels.
1536 if(imulticount<3)
1537 Trigger.SetTopology(1);
1538 else
1539 {
1540 // We should be careful that topologies are sort from
1541 // the less to the more restrictive one.
1542 Trigger.SetTopology(isorttopo[itopocount]);
1543 }
1544 Trigger.ClearFirst();
1545
1546 //
1547 // Start the First Level Trigger simulation
1548 //
1549 if(Lev0!=0)
1550 Lev1=Trigger.FirstLevel();
1551 if(Lev1>0) {
1552 btrigger= 1;
1553 ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++;
1554 }
1555
1556 Lev0=1;
1557 Int_t NumImages = Lev1;
1558 if(Lev1==0 && (Write_All_Images || btrigger)){
1559 btrigger= 1;
1560 NumImages=1;
1561 Lev0=0;
1562 }
1563
1564 for (Int_t ii=0;ii<NumImages;ii++){
1565 if (Write_McTrig){
1566 McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
1567 McTrig[iconcount]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
1568 Trigger.GetMapDiskriminator(trigger_map);
1569 McTrig[iconcount]->SetMapPixels(trigger_map,ii);
1570 }
1571 //
1572 // fill inside the class fadc the member output
1573 //
1574
1575 fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
1576
1577 if( Write_RawEvt ){
1578 //
1579 // Fill the header of this event
1580 //
1581
1582 EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
1583 // fill pixel information
1584 if (Lev1){
1585 for(UInt_t i=0;i<camgeom.GetNumPixels();i++){
1586 for (j=0;j<FADC_SLICES;j++){
1587 fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
1588 }
1589 EvtData[iconcount]->AddPixel(i,fadcValues,0);
1590 }
1591 }
1592 }
1593 }
1594 //
1595 // Increase counter of analised trigger conditions
1596 //
1597 iconcount++;
1598 }
1599 }
1600 else{
1601 break;
1602 }
1603 }
1604 if (!btrigger) break;
1605 }
1606 if (btrigger){
1607
1608 //
1609 // fill the MMcEvt with all information
1610 //
1611
1612 if (Write_McEvt) {
1613 Float_t ftime, ltime;
1614 if (reflector_file_version<6){
1615 mcevth.get_times(&ftime, &ltime);
1616 McEvt->Fill( -1.0,
1617 (UShort_t) mcevth.get_primary() ,
1618 mcevth.get_energy(),
1619 -1.0,
1620 -1.0,
1621 -1.0,
1622 mcevth.get_theta(),
1623 mcevth.get_phi(),
1624 mcevth.get_core(),
1625 mcevth.get_coreX(),
1626 mcevth.get_coreY(),
1627 impactD,
1628 mcevth_2.get_theta_CT(),
1629 mcevth_2.get_phi_CT(),
1630 ftime,
1631 ltime,
1632 0,
1633 0,
1634 0,
1635 0,
1636 0,
1637 0,
1638 0,
1639 (UInt_t)mcevth.get_CORSIKA(),
1640 (UInt_t)mcevth.get_AtmAbs(),
1641 (UInt_t)mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot(),
1642 (UInt_t) ncph,
1643 (UInt_t) inumphe,
1644 (UInt_t) inumphensb+inumphe);
1645 }
1646 else{
1647 Float_t Nmax, t0, tmax, a, b, c, chi2;
1648 mcevth.get_times(&ftime, &ltime);
1649 chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
1650 McEvt->Fill( mcevth_2.get_evt_number(),
1651 (UShort_t) mcevth_2.get_primary() ,
1652 mcevth_2.get_energy(),
1653 mcevth_2.get_thick0(),
1654 mcevth_2.get_first_target(),
1655 mcevth_2.get_z_first_int(),
1656 mcevth_2.get_theta(),
1657 mcevth_2.get_phi(),
1658 mcevth_2.get_core(),
1659 mcevth_2.get_coreX(),
1660 mcevth_2.get_coreY(),
1661 impactD,
1662 mcevth_2.get_theta_CT(),
1663 mcevth_2.get_phi_CT(),
1664 ftime,
1665 ltime,
1666 Nmax,
1667 t0,
1668 tmax,
1669 a,
1670 b,
1671 c,
1672 chi2,
1673 (UInt_t)mcevth_2.get_CORSIKA(),
1674 (UInt_t)mcevth_2.get_AtmAbs(),
1675 (UInt_t)mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot(),
1676 (UInt_t) ncph,
1677 (UInt_t) inumphe,
1678 (UInt_t) inumphensb+inumphe);
1679 }
1680 }
1681 // Fill the Tree with the current leaves of each branch
1682 i=EvtTree.Fill() ;
1683
1684 // Clear the branches
1685 if(Write_McTrig){
1686 for(i=0;i<icontrigger;i++){
1687 McTrig[i]->Clear() ;
1688 }
1689 }
1690 if( Write_RawEvt ){
1691 for(i=0;i<icontrigger;i++){
1692 EvtHeader[i]->Clear() ;
1693 EvtData[i]->DeletePixels();
1694 }
1695 }
1696 if (Write_McEvt)
1697 McEvt->Clear() ;
1698 }
1699 }
1700 // We study a single trigger condition
1701 else {
1702
1703 // Setting trigger conditions
1704 Trigger.SetMultiplicity(Trigger_multiplicity);
1705 Trigger.SetTopology(Trigger_topology);
1706 for (int i=0;i<CAMERA_PIXELS;i++)
1707 fpixelthres[i]=qThreshold[i];
1708 Trigger.SetThreshold(fpixelthres);
1709
1710 Trigger.Diskriminate() ;
1711
1712 //
1713 // look if in all the signals in the trigger signal branch
1714 // is a possible Trigger. Therefore we habe to diskriminate all
1715 // the simulated analog signals (Method Diskriminate in class
1716 // MTrigger). We look simultanously for the moments at which
1717 // there are more than TRIGGER_MULTI pixels above the
1718 // CHANNEL_THRESHOLD.
1719 //
1720
1721 Lev0 = (Short_t) Trigger.ZeroLevel() ;
1722
1723 Lev1 = 0 ;
1724
1725 //
1726 // Start the First Level Trigger simulation
1727 //
1728
1729 if ( Lev0 > 0 || Write_All_Images) {
1730 Lev1= Trigger.FirstLevel();
1731 }
1732 if (Lev1>0){
1733 ++ntrigger;
1734 }
1735 Int_t NumImages = Lev1;
1736 Lev0=1;
1737 if (Lev1==0 && Write_All_Images){
1738 NumImages=1;
1739 Lev0=0;
1740 }
1741
1742 for(Int_t ii=0;ii<NumImages;ii++){
1743 // Loop over different level one triggers
1744
1745 //
1746 // fill inside class fadc the member output
1747 //
1748 fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
1749
1750 if (Write_McTrig){
1751 McTrig[0]->SetFirstLevel ((ii+1)*Lev0);
1752 McTrig[0]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
1753 Trigger.GetMapDiskriminator(trigger_map);
1754 McTrig[0]->SetMapPixels(trigger_map,ii);
1755 }
1756
1757 // Fill Evt information
1758
1759 if (Write_RawEvt){
1760
1761 //
1762 // Fill the header of this event
1763 //
1764
1765 EvtHeader[0]->FillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ;
1766
1767 // fill pixel information
1768
1769 if (Lev1){
1770 for(UInt_t i=0;i<camgeom.GetNumPixels();i++){
1771 for (j=0;j<FADC_SLICES;j++){
1772 fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
1773 }
1774 EvtData[0]->AddPixel(i,fadcValues,0);
1775 }
1776 }
1777 }
1778 //
1779 // fill the MMcEvt with all information
1780 //
1781
1782 if (Write_McEvt){
1783 Float_t ftime, ltime;
1784 if (reflector_file_version<6){
1785 mcevth.get_times(&ftime, &ltime);
1786 McEvt->Fill( -1.0,
1787 (UShort_t) mcevth.get_primary() ,
1788 mcevth.get_energy(),
1789 -1.0,
1790 -1.0,
1791 -1.0,
1792 mcevth.get_theta(),
1793 mcevth.get_phi(),
1794 mcevth.get_core(),
1795 mcevth.get_coreX(),
1796 mcevth.get_coreY(),
1797 impactD,
1798 mcevth_2.get_theta_CT(),
1799 mcevth_2.get_phi_CT(),
1800 ftime,
1801 ltime,
1802 0,
1803 0,
1804 0,
1805 0,
1806 0,
1807 0,
1808 0,
1809 (UInt_t)mcevth.get_CORSIKA(),
1810 (UInt_t)mcevth.get_AtmAbs(),
1811 (UInt_t)mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot(),
1812 (UInt_t) ncph,
1813 (UInt_t) inumphe,
1814 (UInt_t) inumphensb+inumphe);
1815 }
1816 else{
1817 Float_t Nmax, t0, tmax, a, b, c, chi2;
1818 mcevth.get_times(&ftime, &ltime);
1819 chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
1820 McEvt->Fill( mcevth_2.get_evt_number(),
1821 (UShort_t) mcevth_2.get_primary() ,
1822 mcevth_2.get_energy(),
1823 mcevth_2.get_thick0(),
1824 mcevth_2.get_first_target(),
1825 mcevth_2.get_z_first_int(),
1826 mcevth_2.get_theta(),
1827 mcevth_2.get_phi(),
1828 mcevth_2.get_core(),
1829 mcevth_2.get_coreX(),
1830 mcevth_2.get_coreY(),
1831 impactD,
1832 mcevth_2.get_theta_CT(),
1833 mcevth_2.get_phi_CT(),
1834 ftime,
1835 ltime,
1836 Nmax,
1837 t0,
1838 tmax,
1839 a,
1840 b,
1841 c,
1842 chi2,
1843 (UInt_t)mcevth_2.get_CORSIKA(),
1844 (UInt_t)mcevth_2.get_AtmAbs(),
1845 (UInt_t)mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot(),
1846 (UInt_t) ncph,
1847 (UInt_t) inumphe,
1848 (UInt_t) inumphensb+inumphe);
1849 }
1850 }
1851 // We don not count photons out of the camera.
1852
1853
1854 //
1855 // write it out to the file outfile
1856 //
1857
1858 EvtTree.Fill() ;
1859
1860 //
1861 // if a first level trigger occurred, then
1862 // 1. do some other stuff (not implemented)
1863 // 2. start the gui tool
1864
1865 if(FADC_Scan){
1866 if ( Lev0 > 0 ) {
1867 fadc.ShowSignal( McEvt, (Float_t) 60. ) ;
1868 }
1869 }
1870
1871 if(Trigger_Scan){
1872 if ( Lev0 > 0 ) {
1873 Trigger.ShowSignal(McEvt) ;
1874 }
1875 }
1876
1877 // clear all
1878 if (Write_RawEvt) EvtHeader[0]->Clear() ;
1879 if (Write_RawEvt) EvtData[0]->DeletePixels();
1880 if (Write_McEvt) McEvt->Clear() ;
1881 }
1882 if (Write_McTrig) McTrig[0]->Clear() ;
1883 }
1884
1885#ifdef __DEBUG__
1886 printf("\n");
1887
1888 for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1889
1890 for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1891
1892 if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1893
1894 if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1895
1896 printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1897 (int)pixels[ici][icj][PIXNUM],
1898 pixels[ici][icj][PIXX],
1899 pixels[ici][icj][PIXY],
1900 fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1901
1902 }
1903 }
1904 }
1905 }
1906
1907 for (int i=0; i<camgeom.GetNumPixels(); ++i) {
1908 printf("%d (%d): ", i, npixneig[i]);
1909 for (j=0; j<npixneig[i]; ++i)
1910 printf(" %d", pixneig[i][j]);
1911 printf("\n");
1912 }
1913
1914#endif // __DEBUG__
1915
1916
1917 // look for the next event
1918
1919 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1920
1921 } // end while there is a next event
1922
1923 if( !isA( flag, FLAG_END_OF_RUN )){
1924 error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1925 break;
1926 }
1927 else { // found end of run
1928 ntshow += nshow;
1929 log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1930
1931 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1932
1933 if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1934 log(SIGNATURE, "End of file . . .\n");
1935 still_in_loop = FALSE;
1936 log(SIGNATURE, "Reading ASCII files at the end of the reflector file. . .\n");
1937 read_ascii(inputfile);
1938 break;
1939 if ((! Data_From_STDIN) && ( !feof(inputfile) )){
1940
1941 // we have concatenated input files.
1942 // get signature of the next part and check it.
1943
1944 if((reflector_file_version=check_reflector_file( inputfile ))==FALSE){
1945 // exit(1);
1946 log(SIGNATURE, "Next file is not recognised as a reflector file.\n");
1947 log(SIGNATURE, "Stopping ...\n");
1948 break;
1949 }
1950
1951 }
1952
1953 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1954
1955 } // end if found end of file
1956 } // end if found end of run
1957
1958 } // end if else found start of run
1959 } // end big while loop
1960
1961 //<@ Finally we should fill th McRunHeader
1962
1963 Float_t heights[10];
1964 time_t ltime;
1965 Float_t ftime;
1966 Float_t rnum;
1967
1968 get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
1969 get_teles_axis(&telestheta,&telesphi);
1970 get_source_off(&sofftheta,&soffphi);
1971 mcevth.get_theta_range(&shthetamin, &shthetamax);
1972 mcevth.get_phi_range(&shphimin,&shphimax);
1973 corsika=UInt_t(mcevth.get_VersionPGM()*1000);
1974 if(!Trigger_Loop) icontrigger=0;
1975 time (&ltime);
1976 ftime = ((Float_t)ltime)/1000;
1977 for (int i=0; i< 10;i++)
1978 heights[i]=mcevth.get_HeightLev (i);
1979 rnum=mcevth.get_RunNumber();
1980
1981 McRunHeader->Fill(rnum,
1982 (UInt_t) 0,
1983 mcevth.get_DateRun(),
1984 ftime,
1985 icontrigger,
1986 !Write_All_Images,
1987 Write_McEvt,
1988 Write_McTrig,
1989 Write_McFADC,
1990 Write_RawEvt,
1991 addElecNoise,
1992 CAMERA_PIXELS,
1993 (UInt_t)ntshow,
1994 (UInt_t)ntrigger,
1995 0,
1996 sfRaH,
1997 sfRaM,
1998 sfRaS,
1999 sfDeD,
2000 sfDeM,
2001 sfDeS,
2002 meanNSB,
2003 telestheta,
2004 telesphi,
2005 sofftheta,
2006 soffphi,
2007 shthetamax,
2008 shthetamin,
2009 shphimax,
2010 shphimin,
2011 mcevth.get_CWaveLower(),
2012 mcevth.get_CWaveUpper(),
2013 mcevth.get_slope(),
2014 1,
2015 heights,
2016 corsika,
2017 (UInt_t)(REFL_VERSION_A*100),
2018 (UInt_t)(VERSION*100),
2019 0);
2020
2021 // Fill the Header Tree with the current leaves of each branch
2022 HeaderTree.Fill() ;
2023
2024 //++
2025 // put the Event to the root file
2026 //--
2027
2028 outfile_temp.Write() ;
2029 outfile_temp.Close() ;
2030
2031 // close input file
2032
2033 log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
2034 ntshow, ntcph );
2035 datafile<<ntshow<<" event(s), with a total of "<<ntcph<<" C.photons"<<endl;
2036 if (Trigger_Loop){
2037 log( SIGNATURE, "Trigger Mode. \n");
2038 log( SIGNATURE, "Fraction of triggers: \n");
2039 datafile<<"Fraction of triggers: "<<endl;
2040 for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
2041 for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
2042 for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
2043 log( SIGNATURE, "Thres %5.1f, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
2044 fthrescount,imulticount,isorttopo[itopocount],((float)ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0), ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop], ntshow);
2045 datafile<<"Thres "<<fthrescount<<", Multi "<<imulticount<<", Topo"<<isorttopo[itopocount]<<": ";
2046 datafile<<((float)ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0)<<"% ("<<ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]<<" out of "<<ntshow<<")"<<endl;
2047 }
2048 }
2049 }
2050 }
2051 else{
2052 log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
2053 ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
2054 datafile<<"Fraction of triggers: "<<((float)ntrigger) / ((float)ntshow) * 100.0<<" ("<<ntrigger<<" out of "<<ntshow<<" )"<<endl;
2055 }
2056
2057 // close files
2058
2059 log( SIGNATURE, "Closing files\n" );
2060
2061 if( ! Data_From_STDIN ){
2062 fclose( inputfile );
2063 }
2064 datafile.close();
2065
2066 // program finished
2067
2068 log( SIGNATURE, "Done.\n");
2069
2070 return( 0 );
2071}
2072//!@}
2073
2074// @T \newpage
2075
2076//!@subsection Functions definition.
2077
2078//!-----------------------------------------------------------
2079// @name present
2080//
2081// @desc Make some presentation
2082//
2083// @date Sat Jun 27 05:58:56 MET DST 1998
2084//------------------------------------------------------------
2085// @function
2086
2087//!@{
2088void
2089present(void)
2090{
2091 cout << "##################################################\n"
2092 << SIGNATURE << '\n' << '\n'
2093 << "Processor of the reflector output\n"
2094 << "J C Gonzalez, Jun 1998\n"
2095 << "##################################################\n\n"
2096 << flush ;
2097}
2098//!@}
2099
2100
2101//!-----------------------------------------------------------
2102// @name usage
2103//
2104// @desc show help
2105//
2106// @date Tue Dec 15 16:23:30 MET 1998
2107//------------------------------------------------------------
2108// @function
2109
2110//!@{
2111void
2112usage(void)
2113{
2114 present();
2115 cout << "\nusage ::\n\n"
2116 << "\t camera "
2117 << " [ -@ paramfile ] "
2118 << " [ -h ] "
2119 << "\n\n or \n\n"
2120 << "\t camera < paramfile"
2121 << "\n\n";
2122 exit(0);
2123}
2124//!@}
2125
2126
2127//!-----------------------------------------------------------
2128// @name log
2129//
2130// @desc function to send log information
2131//
2132// @var funct Name of the caller function
2133// @var fmt Format to be used (message)
2134// @var ... Other information to be shown
2135//
2136// @date Sat Jun 27 05:58:56 MET DST 1998
2137//------------------------------------------------------------
2138// @function
2139
2140//!@{
2141void
2142log(const char *funct, char *fmt, ...)
2143{
2144 va_list args;
2145
2146 // Display the name of the function that called error
2147 printf("[%s]: ", funct);
2148
2149 // Display the remainder of the message
2150 va_start(args, fmt);
2151 vprintf(fmt, args);
2152 va_end(args);
2153}
2154//!@}
2155
2156
2157//!-----------------------------------------------------------
2158// @name error
2159//
2160// @desc function to send an error message, and abort the program
2161//
2162// @var funct Name of the caller function
2163// @var fmt Format to be used (message)
2164// @var ... Other information to be shown
2165//
2166// @date Sat Jun 27 05:58:56 MET DST 1998
2167//------------------------------------------------------------
2168// @function
2169
2170//!@{
2171void
2172error(const char *funct, char *fmt, ...)
2173{
2174 va_list args;
2175
2176 // Display the name of the function that called error
2177 fprintf(stdout, "ERROR in %s: ", funct);
2178
2179 // Display the remainder of the message
2180 va_start(args, fmt);
2181 vfprintf(stdout, fmt, args);
2182 va_end(args);
2183
2184 perror(funct);
2185
2186 exit(1);
2187}
2188//!@}
2189
2190
2191//!-----------------------------------------------------------
2192// @name isA
2193//
2194// @desc returns TRUE(FALSE), if the flag is(is not) the given
2195//
2196// @var s1 String to be searched
2197// @var flag Flag to compare with string s1
2198// @return TRUE: both strings match; FALSE: oth.
2199//
2200// @date Wed Jul 8 15:25:39 MET DST 1998
2201//------------------------------------------------------------
2202// @function
2203
2204//!@{
2205int
2206isA( char * s1, const char * flag ) {
2207 return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
2208}
2209//!@}
2210
2211
2212//!-----------------------------------------------------------
2213// @name read_ct_file
2214//
2215// @desc read CT definition file
2216//
2217// @date Sat Jun 27 05:58:56 MET DST 1998
2218//------------------------------------------------------------
2219// @function
2220
2221//!@{
2222void
2223read_ct_file(void)
2224{
2225 char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
2226 char token[ITEM_MAX_LENGTH]; //@< a single token
2227 int i, j; //@< dummy counters
2228
2229 log( "read_ct_file", "start.\n" );
2230
2231 ifstream ctin ( ct_filename );
2232
2233 if ( ctin.bad() )
2234 error( "read_ct_file",
2235 "Cannot open CT def. file: %s\n", ct_filename );
2236
2237 // loop till the "end" directive is reached
2238
2239 while (!ctin.eof()) {
2240
2241 // get line from stdin
2242
2243 ctin.getline(line, LINE_MAX_LENGTH);
2244
2245 // look for each item at the beginning of the line
2246
2247 for (i=0; i<=define_mirrors; i++)
2248 if (strstr(line, CT_ITEM_NAMES[i]) == line)
2249 break;
2250
2251 // if it is not a valid line, just ignore it
2252
2253 if (i == define_mirrors+1)
2254 continue;
2255
2256 // case block for each directive
2257
2258 switch ( i ) {
2259
2260 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2261
2262 // get focal distance
2263
2264 sscanf(line, "%s %d", token, &ct_Type);
2265
2266 log( "read_ct_file", "<Type of Telescope>: %s\n",
2267 ((ct_Type==0) ? "CT1" : "MAGIC") );
2268
2269 break;
2270
2271 case focal_distance: // <focal distance> [cm]
2272
2273 // get focal distance
2274
2275 sscanf(line, "%s %f", token, &ct_Focal_mean);
2276
2277 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2278
2279 break;
2280
2281 case focal_std: // s(focal distance) [cm]
2282
2283 // get focal distance
2284
2285 sscanf(line, "%s %f", token, &ct_Focal_std);
2286
2287 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2288
2289 break;
2290
2291 case point_spread: // <point spread> [cm]
2292
2293 // get point spread
2294
2295 sscanf(line, "%s %f", token, &ct_PSpread_mean);
2296
2297 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2298
2299 break;
2300
2301 case point_std: // s(point spread) [cm]
2302
2303 // get point spread
2304
2305 sscanf(line, "%s %f", token, &ct_PSpread_std);
2306
2307 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2308
2309 break;
2310
2311 case adjustment_dev: // s(adjustment_dev) [cm]
2312
2313 // get point spread
2314
2315 sscanf(line, "%s %f", token, &ct_Adjustment_std);
2316
2317 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2318
2319 break;
2320
2321 case black_spot: // radius of the black spot in the center [cm]
2322
2323 // get black spot radius
2324
2325 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2326
2327 log( "read_ct_file", "Radius of the black spots: %f cm\n",
2328 ct_BlackSpot_rad);
2329
2330 break;
2331
2332 case r_mirror: // radius of the mirrors [cm]
2333
2334 // get radius of mirror
2335
2336 sscanf(line, "%s %f", token, &ct_RMirror);
2337
2338 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2339
2340 break;
2341
2342 case n_mirrors: // number of mirrors
2343
2344 // get the name of the output_file from the line
2345
2346 sscanf(line, "%s %d", token, &ct_NMirrors);
2347
2348 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2349
2350 break;
2351
2352 case camera_width: // camera width [cm]
2353
2354 // get the name of the ct_file from the line
2355
2356 sscanf(line, "%s %f", token, &ct_CameraWidth);
2357
2358 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2359
2360 break;
2361
2362 case n_pixels: // number of pixels
2363
2364 // get the name of the output_file from the line
2365
2366 sscanf(line, "%s %d", token, &ct_NPixels);
2367
2368 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2369
2370 break;
2371
2372 case n_centralpixels: // number of central pixels
2373
2374 // get the name of the output_file from the line
2375
2376 sscanf(line, "%s %d", token, &ct_NCentralPixels);
2377
2378 log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
2379
2380 break;
2381
2382 case n_gappixels: // number of gap pixels
2383
2384 // get the name of the output_file from the line
2385
2386 sscanf(line, "%s %d", token, &ct_NGapPixels);
2387
2388 log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
2389
2390 break;
2391
2392 case pixel_width: // pixel width [cm]
2393
2394 // get the name of the ct_file from the line
2395
2396 sscanf(line, "%s %f", token, &ct_PixelWidth);
2397
2398 ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
2399
2400 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2401
2402 break;
2403
2404 case define_mirrors: // read table with the parameters of the mirrors
2405
2406 log( "read_ct_file", "Table of mirrors data:\n" );
2407
2408 // check whether the number of mirrors was already set
2409
2410 if ( ct_NMirrors == 0 )
2411 error( "read_ct_file", "NMirrors was not set.\n" );
2412
2413 // allocate memory for paths list
2414
2415 log( "read_ct_file", "Allocating memory for ct_data\n" );
2416
2417 ct_data = new float*[ct_NMirrors];
2418
2419 for (i=0; i<ct_NMirrors; i++)
2420 ct_data[i] = new float[CT_NDATA];
2421
2422 // read data
2423
2424 log( "read_ct_file", "Reading mirrors data...\n" );
2425
2426 for (i=0; i<ct_NMirrors; i++)
2427 for (j=0; j<CT_NDATA; j++)
2428 ctin >> ct_data[i][j];
2429
2430 break;
2431
2432 } // switch ( i )
2433
2434 } // end while
2435
2436 // end
2437
2438 log( "read_ct_file", "done.\n" );
2439
2440 return;
2441}
2442//!@}
2443
2444
2445//!-----------------------------------------------------------
2446// @name read_QE
2447//
2448// @desc read QE data
2449//
2450// @date thu 5 17:59:57 CEST 2002
2451//------------------------------------------------------------
2452// @function
2453
2454//!@{
2455void
2456read_QE(void){
2457 ifstream qefile;
2458 char line[LINE_MAX_LENGTH];
2459 int i, j, icount;
2460 float qe;
2461
2462 //------------------------------------------------------------
2463 // second, pixels' QE
2464
2465 // try to open the file
2466
2467 log("read_QE", "Opening the file \"%s\" . . .\n", QE_FILE);
2468
2469 qefile.open( QE_FILE );
2470
2471 // if it is wrong or does not exist, exit
2472
2473 if ( qefile.bad() )
2474 error( "read_QE", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2475
2476 // read file
2477
2478 log("read_QE", "Reading data . . .\n");
2479
2480 i=-1;
2481 icount = 0;
2482
2483 while ( ! qefile.eof() ) {
2484
2485 // get line from the file
2486
2487 qefile.getline(line, LINE_MAX_LENGTH);
2488
2489 // skip if comment
2490
2491 if ( *line == '#' )
2492 continue;
2493
2494 // if it is the first valid value, it is the number of QE data points
2495
2496 if ( i < 0 ) {
2497
2498 // get the number of datapoints
2499
2500 sscanf(line, "%d", &pointsQE);
2501
2502 // allocate memory for the table of QEs
2503
2504 QE = new float ** [ct_NPixels];
2505
2506 for ( i=0; i<ct_NPixels; ++i ) {
2507 QE[i] = new float * [2];
2508 QE[i][0] = new float[pointsQE];
2509 QE[i][1] = new float[pointsQE];
2510 }
2511
2512 QElambda = new float [pointsQE];
2513
2514 for ( i=0; i<pointsQE; ++i ) {
2515 qefile.getline(line, LINE_MAX_LENGTH);
2516 sscanf(line, "%f", &QElambda[i]);
2517 }
2518
2519 i=0;
2520
2521 continue;
2522 }
2523
2524 // get the values (num-pixel, num-datapoint, QE-value)
2525
2526 if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2527 break;
2528
2529 if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2530 ((j-1) < pointsQE) && ((j-1) > -1) ) {
2531 QE[i-1][0][j-1] = QElambda[j-1];
2532 QE[i-1][1][j-1] = qe;
2533 }
2534
2535 if ( i > ct_NPixels) break;
2536
2537 icount++;
2538
2539 }
2540
2541 if(icount/pointsQE < ct_NPixels){
2542 error( "read_QE", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2543 icount/pointsQE, ct_NPixels );
2544 }
2545
2546 // close file
2547
2548 qefile.close();
2549
2550 // test QE
2551
2552 for(icount=0; icount< ct_NPixels; icount++){
2553 for(i=0; i<pointsQE; i++){
2554 if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
2555 QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
2556 error( "read_QE", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2557 icount, i, QE[icount][0][i], QE[icount][1][i] );
2558 }
2559 }
2560 }
2561
2562 // end
2563
2564 log("read_QE", "Done.\n");
2565}
2566//!@}
2567
2568//!-----------------------------------------------------------
2569// @name read_WC
2570//
2571// @desc read WC data
2572//
2573// @date thu 5 17:59:57 CEST 2002
2574//------------------------------------------------------------
2575// @function
2576
2577//!@{
2578void
2579read_WC(void){
2580 ifstream wcfile;
2581 char line[LINE_MAX_LENGTH];
2582 int i;
2583
2584 //------------------------------------------------------------
2585 // Read Light Guides data
2586
2587 // try to open the file
2588
2589 log("read_QE", "Opening the file \"%s\" . . .\n", WC_FILE);
2590
2591 wcfile.open( WC_FILE );
2592
2593 // if it is wrong or does not exist, exit
2594
2595 if ( wcfile.bad() )
2596 error( "read_WC", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2597
2598 // read file
2599
2600 log("read_WC", "Reading data . . .\n");
2601
2602 // get line from the file
2603
2604 wcfile.getline(line, LINE_MAX_LENGTH);
2605
2606 // get the number of datapoints
2607
2608 sscanf(line, "%d", &pointsWC);
2609
2610 // allocate memory for the table of QEs
2611
2612 WC = new float * [2];
2613 WC[0] = new float[pointsWC];
2614 WC[1] = new float[pointsWC];
2615
2616 for ( i=0; i<pointsWC; ++i ) {
2617 wcfile.getline(line, LINE_MAX_LENGTH);
2618 sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
2619 }
2620
2621 // close file
2622
2623 wcfile.close();
2624
2625 // read
2626
2627 log("read_WC", "Done.\n");
2628}
2629//!@}
2630
2631//!-----------------------------------------------------------
2632// @name read_pixels
2633//
2634// @desc read pixels data
2635//
2636// @date Fri Mar 12 16:33:34 MET 1999
2637//------------------------------------------------------------
2638// @function
2639
2640//!@{
2641void
2642read_pixels(struct camera *pcam)
2643{
2644 ifstream qefile;
2645 char line[LINE_MAX_LENGTH];
2646 int n, i, j, icount;
2647 float qe;
2648
2649 //------------------------------------------------------------
2650 // first, pixels' coordinates
2651
2652 pcam->inumpixels = ct_NPixels;
2653 pcam->inumcentralpixels = ct_NCentralPixels;
2654 pcam->inumgappixels = ct_NGapPixels;
2655 pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
2656 pcam->dpixdiameter_cm = ct_PixelWidth;
2657
2658 // initialize pixel numbers
2659
2660 pixary = new float* [2*ct_NCentralPixels];
2661 for ( i=0; i<2*ct_NCentralPixels; ++i )
2662 pixary[i] = new float[2];
2663
2664 pixneig = new int* [ct_NCentralPixels];
2665 for ( i=0; i<ct_NCentralPixels; ++i ) {
2666 pixneig[i] = new int[6];
2667 for ( j=0; j<6; ++j )
2668 pixneig[i][j] = -1;
2669 }
2670
2671 npixneig = new int[ct_NCentralPixels];
2672 for ( i=0; i<ct_NCentralPixels; ++i )
2673 npixneig[i] = 0;
2674
2675 // generate all coordinates
2676
2677 igen_pixel_coordinates(pcam);
2678
2679
2680 // calculate tables of neighbours
2681
2682#ifdef __DEBUG__
2683 for ( n=0 ; n<ct_NPixels ; ++n ) {
2684 cout << "Para el pixel " << n << ": ";
2685 for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2686 if ( pixels_are_neig(n,i) == TRUE ) {
2687 pixneig[n][npixneig[n]] = i;
2688 pixneig[i][npixneig[i]] = n;
2689 cout << i << ' ';
2690 ++npixneig[n];
2691 ++npixneig[i];
2692 }
2693 }
2694 cout << endl << flush;
2695 }
2696#else // ! __DEBUG__
2697 for ( n=0 ; n<ct_NCentralPixels ; ++n )
2698 for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2699 if ( pixels_are_neig(n,i) == TRUE ) {
2700 pixneig[n][npixneig[n]] = i;
2701 pixneig[i][npixneig[i]] = n;
2702 ++npixneig[n];
2703 ++npixneig[i];
2704 }
2705#endif // ! __DEBUG__
2706
2707#ifdef __DEBUG__
2708 for ( n=0 ; n<ct_NPixels ; ++n ) {
2709 cout << n << ':';
2710 for ( j=0; j<npixneig[n]; ++j)
2711 cout << ' ' << pixneig[n][j];
2712 cout << endl << flush;
2713 }
2714#endif // __DEBUG__
2715
2716 //------------------------------------------------------------
2717 // second, pixels' QE
2718
2719 // try to open the file
2720
2721 log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
2722
2723 qefile.open( QE_FILE );
2724
2725 // if it is wrong or does not exist, exit
2726
2727 if ( qefile.bad() )
2728 error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2729
2730 // read file
2731
2732 log("read_pixels", "Reading data . . .\n");
2733
2734 i=-1;
2735 icount = 0;
2736
2737 while ( ! qefile.eof() ) {
2738
2739 // get line from the file
2740
2741 qefile.getline(line, LINE_MAX_LENGTH);
2742
2743 // skip if comment
2744
2745 if ( *line == '#' )
2746 continue;
2747
2748 // if it is the first valid value, it is the number of QE data points
2749
2750 if ( i < 0 ) {
2751
2752 // get the number of datapoints
2753
2754 sscanf(line, "%d", &pointsQE);
2755
2756 // allocate memory for the table of QEs
2757
2758 QE = new float ** [ct_NPixels];
2759
2760 for ( i=0; i<ct_NPixels; ++i ) {
2761 QE[i] = new float * [2];
2762 QE[i][0] = new float[pointsQE];
2763 QE[i][1] = new float[pointsQE];
2764 }
2765
2766 QElambda = new float [pointsQE];
2767
2768 for ( i=0; i<pointsQE; ++i ) {
2769 qefile.getline(line, LINE_MAX_LENGTH);
2770 sscanf(line, "%f", &QElambda[i]);
2771 }
2772
2773 i=0;
2774
2775 continue;
2776 }
2777
2778 // get the values (num-pixel, num-datapoint, QE-value)
2779
2780 if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2781 break;
2782
2783 if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2784 ((j-1) < pointsQE) && ((j-1) > -1) ) {
2785 QE[i-1][0][j-1] = QElambda[j-1];
2786 QE[i-1][1][j-1] = qe;
2787 }
2788
2789 if ( i > ct_NPixels) break;
2790
2791 icount++;
2792
2793 }
2794
2795 if(icount/pointsQE < ct_NPixels){
2796 error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2797 icount/pointsQE, ct_NPixels );
2798 }
2799
2800 // close file
2801
2802 qefile.close();
2803
2804 // test QE
2805
2806 for(icount=0; icount< ct_NPixels; icount++){
2807 for(i=0; i<pointsQE; i++){
2808 if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
2809 QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
2810 error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2811 icount, i, QE[icount][0][i], QE[icount][1][i] );
2812 }
2813 }
2814 }
2815
2816 // end
2817
2818 log("read_pixels", "Done.\n");
2819
2820}
2821//!@}
2822
2823//!-----------------------------------------------------------
2824// @name read_ascii
2825//
2826// @desc read ascii configuration files used by the reflector
2827//
2828// @date tue dec 10 17:14:10 CET 2002
2829//------------------------------------------------------------
2830// @function
2831
2832//!@{
2833void
2834read_ascii(FILE *sp)
2835{
2836 cout<<"To be implemented..."<<endl;
2837}
2838
2839//!-----------------------------------------------------------
2840// @name pixels_are_neig
2841//
2842// @desc check whether two pixels are neighbours
2843//
2844// @var pix1 Number of the first pixel
2845// @var pix2 Number of the second pixel
2846// @return TRUE: both pixels are neighbours; FALSE: oth.
2847//
2848// @date Wed Sep 9 17:58:37 MET DST 1998
2849//------------------------------------------------------------
2850// @function
2851
2852//!@{
2853int
2854pixels_are_neig(int pix1, int pix2)
2855{
2856 if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2857 SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2858 > ct_PixelWidth_corner_2_corner )
2859 return ( FALSE );
2860 else
2861 return ( TRUE );
2862}
2863//!@}
2864
2865//!-----------------------------------------------------------
2866// @name igen_pixel_coordinates
2867//
2868// @desc generate the pixel center coordinates
2869//
2870// @var *pcam structure camera containing all the
2871// camera information
2872// @return total number of pixels
2873//
2874// DP
2875//
2876// @date Thu Oct 14 10:41:03 CEST 1999
2877//------------------------------------------------------------
2878// @function
2879
2880//!@{
2881/******** igen_pixel_coordinates() *********************************/
2882
2883int igen_pixel_coordinates(struct camera *pcam) {
2884 /* generate pixel coordinates, return value is number of pixels */
2885
2886 int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2887 float fsegment_fract;
2888 double dtsize;
2889 double dhsize;
2890 double dpsize;
2891 double dxfirst_pix;
2892 double dyfirst_pix;
2893 double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5;
2894 double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5;
2895
2896
2897 double dstartx, dstarty; /* for the gap pixels and outer pixels */
2898 int j, nrow;
2899
2900 dpsize = pcam->dpixdiameter_cm;
2901 dtsize = dpsize * sqrt(3.) / 2.;
2902 dhsize = dpsize / 2.;
2903
2904 /* Loop over central pixels to generate co-ordinates */
2905
2906 for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2907
2908 /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2909
2910 pcam->dpixsizefactor[ipixno-1] = 1.;
2911
2912 in = 0;
2913
2914 i = 0;
2915 itot_inside_ring = 0;
2916 iring_no = 0;
2917
2918 /* Calculate the number of pixels out to and including the ring containing pixel number */
2919 /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2920
2921 while (itot_inside_ring == 0){
2922
2923 iN = 3*(i*(i+1)) + 1;
2924
2925 if (ipixno <= iN){
2926 iring_no = i;
2927 itot_inside_ring = iN;
2928 }
2929
2930 i++;
2931 }
2932
2933
2934 /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2935
2936 ipix_in_ring = 0;
2937 for (i = 0; i < iring_no; ++i){
2938
2939 ipix_in_ring = ipix_in_ring + 6;
2940 }
2941
2942 /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2943 /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2944 /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2945
2946 isegment = 0;
2947 fsegment_fract = 0.;
2948 if (iring_no > 0) {
2949
2950 isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2951
2952 fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2953
2954 }
2955
2956 /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2957 /* pixel width (flat to flat) dpsize. */
2958
2959 dxfirst_pix = dpsize*iring_no;
2960 dyfirst_pix = 0.;
2961
2962 /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2963
2964 ddxseg1 = - dhsize*iring_no;
2965 ddyseg1 = dtsize*iring_no;
2966 ddxseg2 = -dpsize*iring_no;
2967 ddyseg2 = 0.;
2968 ddxseg3 = ddxseg1;
2969 ddyseg3 = -ddyseg1;
2970 ddxseg4 = -ddxseg1;
2971 ddyseg4 = -ddyseg1;
2972 ddxseg5 = -ddxseg2;
2973 ddyseg5 = 0.;
2974
2975 /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2976 /* anti-clockwise around the ring by adding the segment to segment vectors. */
2977
2978 switch (isegment) {
2979
2980 case 0:
2981
2982 pcam->dxc[ipixno-1] = 0.;
2983 pcam->dyc[ipixno-1] = 0.;
2984
2985 case 1:
2986 pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2987 pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2988
2989 break;
2990
2991 case 2:
2992
2993 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2994 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2995
2996 break;
2997
2998 case 3:
2999
3000 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
3001 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
3002
3003 break;
3004
3005 case 4:
3006
3007 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
3008 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
3009
3010 break;
3011
3012 case 5:
3013
3014 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
3015 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
3016
3017 break;
3018
3019 case 6:
3020
3021 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
3022 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
3023
3024 break;
3025
3026 default:
3027
3028 fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
3029 return(0);
3030
3031 } /* end switch */
3032
3033 } /* end for */
3034
3035 dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
3036 dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
3037
3038 if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
3039
3040 j = pcam->inumcentralpixels;
3041
3042 for(i=0; i<pcam->inumgappixels; i=i+6){
3043 pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
3044 pcam->dyc[j + i ] = dstarty;
3045 pcam->dpixsizefactor[j + i] = 1.;
3046 pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
3047 pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
3048 pcam->dpixsizefactor[j + i + 1] = 1.;
3049 pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
3050 pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
3051 pcam->dpixsizefactor[j + i+ 2] = 1.;
3052 pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
3053 pcam->dyc[j + i + 3] = dstarty;
3054 pcam->dpixsizefactor[j + i+ 3] = 1.;
3055 pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
3056 pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
3057 pcam->dpixsizefactor[j + i+ 4] = 1.;
3058 pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
3059 pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
3060 pcam->dpixsizefactor[j + i + 5] = 1.;
3061 } /* end for */
3062 } /* end if */
3063
3064 /* generate positions of the outer pixels */
3065
3066 if( pcam->inumbigpixels > 0 ){
3067
3068 j = pcam->inumcentralpixels + pcam->inumgappixels;
3069
3070 for(i=0; i<pcam->inumbigpixels; i++){
3071 pcam->dpixsizefactor[j + i] = 2.;
3072 }
3073
3074 in = 0;
3075
3076 nrow = (int) ceil(dstartx / 2. / dpsize);
3077
3078 while(in < pcam->inumbigpixels){
3079
3080 pcam->dxc[j + in] = dstartx + dpsize;
3081 pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
3082 pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
3083 pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
3084 pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
3085 pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
3086 pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
3087 pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
3088 pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
3089 pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
3090 pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
3091 pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
3092 for(i=1; i<nrow; i++){
3093 pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
3094 pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
3095 pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
3096 pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
3097 pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
3098 pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
3099 pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
3100 pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
3101 pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
3102 pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
3103 pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
3104 pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
3105 }
3106 in = in + 6 * nrow;
3107 dstartx = dstartx + 2. * dpsize;
3108 nrow = nrow + 1;
3109 } /* end while */
3110
3111 } /* end if */
3112
3113 return(pcam->inumpixels);
3114
3115}
3116//!@}
3117
3118//!-----------------------------------------------------------
3119// @name bpoint_is_in_pix
3120//
3121// @desc check if a point (x,y) in camera coordinates is inside a given pixel
3122//
3123// @var *pcam structure camera containing all the
3124// camera information
3125// @var dx, dy point coordinates in centimeters
3126// @var ipixnum pixel number (starting at 0)
3127// @return TRUE if the point is inside the pixel, FALSE otherwise
3128//
3129// DP
3130//
3131// @date Thu Oct 14 16:59:04 CEST 1999
3132//------------------------------------------------------------
3133// @function
3134
3135//!@{
3136
3137/******** bpoint_is_in_pix() ***************************************/
3138
3139#define sqrt13 0.577350269 // = 1./sqrt(3.)
3140#define sqrt3 1.732050807 // = sqrt(3.)
3141
3142int bpoint_is_in_pix(double dx, double dy, MGeomCamMagic *pgeomcam)
3143{
3144 /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
3145 /* the pixel is assumed to be a "closed set" */
3146
3147 /*
3148 a = length of one of the edges of one pixel,
3149 b = half the width of one pixel
3150 */
3151
3152 const int numN = pgeomcam->GetNumPixels();
3153
3154 for (int i=0; i<numN; i++)
3155 {
3156 MGeomPix &pixel = (*pgeomcam)[i];
3157 const double b = pixel.GetR()/2;
3158 const double a = pixel.GetR()/sqrt3;
3159
3160 const double xx = dx - pixel.GetX();
3161 const double yy = dy - pixel.GetY();
3162
3163 if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
3164 ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){
3165
3166 return i; // inside i
3167 }
3168
3169 // return -1; // outside
3170 }
3171
3172 return -1; // outside
3173}
3174
3175//!@}
3176
3177//------------------------------------------------------------
3178// @name dist_r_P
3179//
3180// @desc distance straight line r - point P
3181//
3182// @date Sat Jun 27 05:58:56 MET DST 1998
3183// @function @code
3184//------------------------------------------------------------
3185// dist_r_P
3186//
3187// distance straight line r - point P
3188//------------------------------------------------------------
3189
3190float
3191dist_r_P(float a, float b, float c,
3192 float l, float m, float n,
3193 float x, float y, float z)
3194{
3195 return (
3196 sqrt((SQR((a-x)*m-(b-y)*l) +
3197 SQR((b-y)*n-(c-z)*m) +
3198 SQR((c-z)*l-(a-x)*n))/
3199 (SQR(l)+SQR(m)+SQR(n))
3200 )
3201 );
3202}
3203
3204//------------------------------------------------------------
3205// @name check_reflector_file
3206//
3207// @desc check if a given reflector file has the right signature
3208// @desc return TRUE or FALSE
3209//
3210// @date Mon Feb 14 16:44:21 CET 2000
3211// @function @code
3212//------------------------------------------------------------
3213
3214int check_reflector_file(FILE *infile){
3215
3216 char sign[20]; // auxiliary variable
3217
3218 strcpy(sign, REFL_SIGNATURE_B);
3219
3220 fread( (char *)sign, strlen(REFL_SIGNATURE_B), 1, infile);
3221 if (strcmp(sign, REFL_SIGNATURE_A) == 0){
3222 fread( (char *)sign, 1, 1, infile);
3223 return 4;
3224 }
3225 else if (strcmp(sign, REFL_SIGNATURE_B) == 0){
3226 fread( (char *)sign, 1, 1, infile);
3227 return 5;
3228 }
3229 else if (strcmp(sign, REFL_SIGNATURE_C) == 0){
3230 fread( (char *)sign, 1, 1, infile);
3231 return 6;
3232 }
3233 else {
3234 cout << "ERROR: Signature of .rfl file is not correct\n";
3235 cout << '"' << sign << '"' << '\n';
3236 cout << "should be: " << REFL_SIGNATURE_A <<" or "<< REFL_SIGNATURE_B <<" or " << REFL_SIGNATURE_C <<" or "<< '\n';
3237 return(FALSE);
3238 }
3239
3240
3241}
3242
3243//------------------------------------------------------------
3244// @name lin_interpol
3245//
3246// @desc interpolate linearly between two points returning the
3247// @desc y-value of the result
3248//
3249// @date Thu Feb 17 11:31:32 CET 2000
3250// @function @code
3251//------------------------------------------------------------
3252
3253float lin_interpol( float x1, float y1, float x2, float y2, float x){
3254
3255 if( (x2 - x1)<=0. ){ // avoid division by zero, return average
3256 cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
3257 return((y1+y2)/2.);
3258 }
3259 else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
3260 return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
3261 }
3262}
3263
3264
3265//------------------------------------------------------------
3266// @name produce_phes
3267//
3268// @desc read the photons of an event, pixelize them and simulate QE
3269// @desc return various statistics and the array of Photoelectrons
3270//
3271// @date Mon Feb 14 16:44:21 CET 2000
3272// @function @code
3273//------------------------------------------------------------
3274
3275int produce_phes( FILE *sp, // the input file
3276 MGeomCamMagic *camgeom, // the camera layout
3277 float minwl_nm, // the minimum accepted wavelength
3278 float maxwl_nm, // the maximum accepted wavelength
3279 class MTrigger *trigger, // the generated phes
3280 class MFadc *fadc,
3281 int *itotnphe, // total number of produced photoelectrons
3282 float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
3283 int *incph, // total number of cph read
3284 float *tmin_ns, // minimum arrival time of all phes
3285 float *tmax_ns, // maximum arrival time of all phes
3286 int star // star=0 starfield rotation star=-1 not rotation
3287 ){
3288
3289 static uint i;
3290 static int k, ipixnum;
3291 static float cx, cy, wl, qe, t;
3292 static float cu, cv, cw;
3293 static MCCphoton photon;
3294 static float **qept;
3295 static char flag[SIZE_OF_FLAGS + 1];
3296 static float radius_mm;
3297 static float SFR , C3 , C2 , C1;
3298 const double pi = 3.14159;
3299
3300 // reset variables
3301
3302 for ( i=0; i<camgeom->GetNumPixels(); ++i ){
3303
3304 nphe[i] = 0.0;
3305
3306 }
3307
3308 *incph = 0;
3309
3310 radius_mm = camgeom->GetMaxRadius();
3311
3312
3313 //- - - - - - - - - - - - - - - - - - - - - - - - -
3314 // read photons and "map" them into the pixels
3315 //--------------------------------------------------
3316
3317 // initialize CPhoton
3318
3319 photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
3320
3321 // read the photons data
3322
3323
3324 // loop over the photons
3325
3326 while (1) {
3327
3328 fread ( flag, SIZE_OF_FLAGS, 1, sp );
3329
3330 if (isA( flag, FLAG_END_OF_EVENT ))
3331 break;
3332
3333 memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
3334
3335 fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
3336
3337
3338 // increase number of photons
3339
3340 (*incph)++;
3341
3342 // Chceck if photon is inside trigger time range
3343
3344 t = photon.get_t() ;
3345
3346 if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
3347 continue;
3348 }
3349
3350 //
3351 // Pixelization
3352 //
3353
3354 // Where the photon provide(0 from starfield, -1 from Mmcs ???
3355
3356 if ( star == 0 && Starfield_rotate == TRUE )
3357
3358 {
3359 get_teles_axis(&zenith, &azimutal);
3360 // Introduction SFR angle
3361
3362 zenith = (pi/180) * zenith;
3363
3364 azimutal =( pi/180) * azimutal;
3365
3366 C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
3367
3368 C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
3369
3370 C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
3371
3372 SFR = acos( C1/C2 );
3373
3374 if ( sin(azimutal) < 0)
3375 {
3376 SFR = 2 * pi - SFR;
3377 }
3378 else
3379 {
3380 SFR = SFR;
3381 }
3382
3383//final of introduction SFR angle//
3384
3385 cx = photon.get_x() * cos(SFR) - photon.get_y() * sin(SFR);
3386 cy = photon.get_y() * cos(SFR) + photon.get_x() * sin(SFR);
3387
3388 }
3389 else
3390 {
3391 cx = photon.get_x();
3392 cy = photon.get_y();
3393
3394 }
3395
3396
3397 // get wavelength
3398
3399 wl = photon.get_wl();
3400
3401 // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n";
3402
3403 // check if photon has valid wavelength and is inside outermost camera radius
3404
3405 if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy)*10 > radius_mm ) ){
3406 continue;
3407
3408 }
3409
3410 ipixnum = bpoint_is_in_pix(cx*10, cy*10, camgeom);
3411
3412 // -1 = the photon is in none of the pixels
3413 // 0 = the phton is in the central pixel, which is not used for trigger
3414 if (ipixnum==-1 || ipixnum==0) {
3415 continue;
3416 }
3417
3418 //+++
3419 // QE simulation
3420 //---
3421
3422 // set pointer to the QE table of the relevant pixel
3423
3424 qept = (float **)QE[ipixnum];
3425
3426 // check if wl is inside table; outside the table, QE is assumed to be zero
3427
3428 if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
3429 continue;
3430
3431 }
3432
3433 // find data point in the QE table (-> k)
3434
3435 k = 1; // start at 1 because the condition was already tested for 0
3436 while (k < pointsQE-1 && qept[0][k] < wl){
3437 k++;
3438 }
3439
3440 // calculate the qe value between 0. and 1.
3441
3442 qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
3443
3444 // Apply incient angular correction due to Light Guides
3445
3446 cu=photon.get_u();
3447 cv=photon.get_v();
3448 cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*pi/180.0;
3449
3450 // find data point in the QE table (-> k)
3451
3452 k = 0;
3453 while (k < pointsWC-1 && WC[0][k] < cw){
3454 k++;
3455 }
3456
3457 // correct the qe with WC data.
3458
3459 qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
3460
3461 // if random > quantum efficiency, reject it
3462
3463 if ( (RandomNumber) > qe ) {
3464 continue;
3465 }
3466
3467 //+++
3468 // The photon has produced a photo electron
3469 //---
3470
3471 // cout << " accepted\n";
3472
3473 // increment the number of photoelectrons in the relevant pixel
3474
3475 nphe[ipixnum] += 1.0;
3476
3477 // store the new photoelectron
3478
3479 fadc->Fill(ipixnum,(t-*tmin_ns) , trigger->FillShow(ipixnum,t-*tmin_ns));
3480
3481 *itotnphe += 1;
3482 }
3483
3484 return(0);
3485
3486}
3487
3488
3489//------------------------------------------------------------
3490// @name produce_nsbrates
3491//
3492// @desc read the starfield file, call produce_phes on it in,
3493// @desc different wavebands, calculate the nsbrates
3494//
3495// @date Mon Feb 14 16:44:21 CET 2000
3496// @function @code
3497//------------------------------------------------------------
3498
3499int produce_nsbrates( char *iname, // the starfield input file name
3500 MGeomCamMagic *camgeom, // camera layout
3501 float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
3502 // the NSB rates in phe/ns for each pixel
3503 ){
3504
3505 uint i, j; // counters
3506 int k, ii; // counters
3507
3508 MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
3509 MFadc flashadc;
3510
3511 static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
3512 WAVEBANDBOUND2,
3513 WAVEBANDBOUND3,
3514 WAVEBANDBOUND4,
3515 WAVEBANDBOUND5,
3516 WAVEBANDBOUND6 };
3517
3518 FILE *infile; // the input file
3519 fpos_t fileposition; // marker on the input file
3520 static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
3521 static MCEventHeader evth; // the event header
3522 static MCEventHeader evth_2; // the event header
3523 static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
3524 int reflector_file_version;
3525 int itnphe; // total number of produced photoelectrons
3526 int itotnphe; // total number of produced photoelectrons after averaging
3527 int incph; // total number of cph read
3528 float tmin_ns; // minimum arrival time of all phes
3529 float tmax_ns; // maximum arrival time of all phes
3530 float integtime_ns; // integration time from the starfield generator
3531
3532 // open input file
3533
3534 log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
3535
3536 infile = fopen( iname, "r" );
3537
3538 // check if the satrfield input file exists
3539
3540 if ( infile == NULL ) {
3541
3542 log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
3543
3544 log( SIGNATURE, "There is not NSB from the Stars\n");
3545
3546 return (0);
3547 }
3548
3549 // get signature, and check it
3550
3551 if((reflector_file_version=check_reflector_file( infile ))==FALSE){
3552 exit(1);
3553 }
3554
3555 // initialize flag
3556
3557 strcpy( cflag, " \0" );
3558
3559 // get flag
3560
3561 fread( cflag, SIZE_OF_FLAGS, 1, infile );
3562
3563 if ( ! feof(infile)){
3564
3565 // reading .rfl file
3566
3567 if(!isA( cflag, FLAG_START_OF_RUN )){
3568 error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
3569 }
3570 else { // found start of run
3571
3572 fread( cflag, SIZE_OF_FLAGS, 1, infile );
3573
3574 if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
3575
3576 // get MCEventHeader
3577
3578 if (reflector_file_version<6)
3579 fread( (char*)&evth, evth.mysize(), 1, infile );
3580 else
3581 fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
3582
3583 if (reflector_file_version<6)
3584 integtime_ns = evth.get_energy();
3585 else
3586 integtime_ns = evth_2.get_energy();
3587
3588 // memorize where we are in the file
3589
3590 if (fgetpos( infile, &fileposition ) != 0){
3591 error( SIGNATURE, "Cannot position in file ...\n");
3592 }
3593
3594 // loop over the wavebands
3595
3596 for(i=0; i<iNUMWAVEBANDS; i++){
3597
3598 // initialize the rate array
3599
3600 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
3601 rate_phepns[j][i] = 0.;
3602 }
3603
3604 itotnphe = 0;
3605
3606 // read the photons and produce the photoelectrons
3607 // - in order to average over the QE simulation, call the
3608 // production function iNUMNSBPRODCALLS times
3609
3610 for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
3611
3612 // position the file pointer to the beginning of the photons
3613
3614 fsetpos( infile, &fileposition);
3615
3616 // produce photoelectrons
3617
3618 k = produce_phes( infile,
3619 camgeom,
3620 wl_nm[i],
3621 wl_nm[i+1],
3622 &trigger, // this is a dummy here
3623 &flashadc, // this is a dummy here
3624 &itnphe,
3625 nphe, // we want this!
3626 &incph,
3627 &tmin_ns,
3628 &tmax_ns,
3629 0 ); // photons from starfield
3630
3631 if( k != 0 ){ // non-zero returnvalue means error
3632 cout << "Exiting.\n";
3633 exit(1);
3634 }
3635
3636 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
3637 rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
3638 }
3639
3640 itotnphe += itnphe;
3641
3642 } // end for(ii=0 ...
3643
3644 fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
3645 wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
3646
3647 } // end for(i=0 ...
3648
3649 }
3650 else{
3651 cout << "Starfield file contains no event.\nExiting.\n";
3652 exit(1);
3653 } // end if( isA ... event
3654 } // end if ( isA ... run
3655 }
3656 else{
3657 cout << "Starfield file contains no run.\nExiting.\n";
3658 exit(1);
3659 }
3660
3661 fclose( infile );
3662
3663 return(0);
3664
3665}
3666
3667
3668//------------------------------------------------------------
3669// @name produce_nsb_phes
3670//
3671// @desc produce the photoelectrons from the nsbrates
3672//
3673// @date Thu Feb 17 17:10:40 CET 2000
3674// @function @code
3675//------------------------------------------------------------
3676
3677int produce_nsb_phes( float *atmin_ns,
3678 float *atmax_ns,
3679 float theta_rad,
3680 struct camera *cam,
3681 float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
3682 float difnsbr_phepns[iMAXNUMPIX],
3683 float extinction[iNUMWAVEBANDS],
3684 float fnpx[iMAXNUMPIX],
3685 MTrigger *trigger,
3686 MFadc *fadc,
3687 int *inphe,
3688 float base_mv[iMAXNUMPIX]){
3689
3690 float simtime_ns; // NSB simulation time
3691 int i, j, k, ii;
3692 float zenfactor; // correction factor calculated from the extinction
3693 int inumnsbphe; // number of photoelectrons caused by NSB
3694
3695 float t;
3696 TRandom random; // Random numbers generator
3697
3698 random.SetSeed ((UInt_t) (RandomNumber*1000));
3699
3700 ii = *inphe; // avoid dereferencing
3701
3702 // check if the arrival times are set; if not generate them
3703
3704 if(*atmin_ns <SIMTIMEOFFSET_NS || *atmin_ns > *atmax_ns){
3705 *atmin_ns = 0.;
3706 *atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
3707
3708 }
3709 else{ // extend the event time window by the given offsets
3710
3711 *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
3712 *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
3713
3714 simtime_ns = *atmax_ns - *atmin_ns;
3715
3716 // make sure the simulated time is long enough for the FADC
3717 // simulation and not too long
3718
3719 if(simtime_ns< TOTAL_TRIGGER_TIME){
3720 *atmin_ns = *atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
3721 *atmax_ns = *atmin_ns + TOTAL_TRIGGER_TIME;
3722 simtime_ns = TOTAL_TRIGGER_TIME;
3723 }
3724
3725 if(simtime_ns> TOTAL_TRIGGER_TIME){
3726 *atmax_ns =*atmin_ns + TOTAL_TRIGGER_TIME;
3727 simtime_ns = TOTAL_TRIGGER_TIME;
3728 }
3729
3730 }
3731
3732 // initialize baselines
3733
3734 for(i=0; i<cam->inumpixels; i++){
3735 base_mv[i] = 0.;
3736 }
3737
3738 // calculate baselines and generate phes
3739
3740 for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
3741
3742 // calculate the effect of the atmospheric extinction
3743
3744 zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
3745
3746 for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
3747
3748 inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
3749 * simtime_ns );
3750
3751 base_mv[j] += inumnsbphe;
3752
3753 // randomize
3754
3755 if (inumnsbphe>0.0){
3756 inumnsbphe = random.Poisson (inumnsbphe );
3757 }
3758
3759 // create the photoelectrons
3760
3761 for(k=0; k < inumnsbphe; k++){
3762
3763 t=(RandomNumber * simtime_ns);
3764
3765 (*fadc).Fill(j,t ,(*trigger).FillNSB(j,t));
3766
3767 ii++; // increment total number of photoelectons
3768
3769 fnpx[j] += 1.; // increment number of photoelectrons in this pixel
3770
3771 }
3772
3773 } // end for(j=0 ...
3774 } // end for(i=0 ...
3775
3776 // finish baseline calculation
3777
3778 for(i=0; i<cam->inumpixels; i++){
3779 base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
3780 }
3781
3782 *inphe = ii; // update the pointer
3783
3784 return(0);
3785}
3786
3787
3788// @endcode
3789
3790
3791//=------------------------------------------------------------
3792//!@subsection Log of this file.
3793
3794//!@{
3795//
3796// $Log: not supported by cvs2svn $
3797// Revision 1.46 2002/11/08 17:51:00 blanch
3798// *** empty log message ***
3799//
3800// Revision 1.45 2002/10/29 17:15:27 blanch
3801// Reading several reflector versions.
3802//
3803// Revision 1.44 2002/10/18 16:53:03 blanch
3804// Modification to read several reflector files.
3805//
3806// Revision 1.43 2002/09/13 10:53:39 blanch
3807// Minor change to remove some undisired comments.
3808//
3809// Revision 1.42 2002/09/09 16:00:49 blanch
3810// Statement has been included to avoid writting to disk MParContainer and MArray.
3811// It has also been added the effect of the WC, the actual values must be added,
3812// once they are measured.
3813//
3814// Revision 1.41 2002/09/04 09:57:42 blanch
3815// Modifications done to use MGeomCam from MARS.
3816//
3817// Revision 1.40 2002/07/16 16:15:22 blanch
3818// A first implementation of the Star field rotation has been done.
3819//
3820// Revision 1.39 2002/04/27 10:48:39 blanch
3821// Some unused varibles have been removed.
3822//
3823// Revision 1.38 2002/03/18 18:44:29 blanch
3824// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
3825//
3826// Revision 1.37 2002/03/18 16:42:20 blanch
3827// The data member fProductionSite of the MMcRunHeader has been set to 0,
3828// which means that the prodution site is unknown.
3829//
3830// Revision 1.35 2002/03/15 15:15:52 blanch
3831// Several modifications to simulate the actual trigger zone.
3832//
3833// Revision 1.34 2002/03/13 18:13:56 blanch
3834// Some changes to fill correctly the new format of MMcRunHeader.
3835//
3836// Revision 1.33 2002/03/04 17:21:48 blanch
3837// Small and not important changes.
3838//
3839// Revision 1.32 2002/02/28 15:04:52 blanch
3840// A small back has been solved. Before, while not using the option
3841// writte_all_images, not all triggered showers were stored. Now it is solved.
3842// For that it is importatn taht the less restrictive trigger option is
3843// checked first.
3844// A new facility has been introduced and now one can choose the step size in
3845// trigger loop mode for the discriminator threshold.
3846// The close-compact topology for less than 3 pixels does not make sense. Before
3847// the program was ignoring that, now it switch to simple neighbour condition.
3848//
3849// Revision 1.31 2002/01/18 17:41:02 blanch
3850// The option of adding noise to all pixels or to not adding the noise
3851// has been added.
3852// We removed the pixels larger than 577. When there were more than one
3853// trigger in one shower, the pixel number was increasing. Now it is
3854// flagged by the variable MMcTrig::fFirstLvlTrig.
3855//
3856// Revision 1.30 2001/11/27 09:49:54 blanch
3857// Fixing bug which was treating wrongly the extension of star photons.
3858//
3859// Revision 1.29 2001/11/14 17:38:23 blanch
3860// Sveral changes have been done:
3861// - bpoint_in_in_pixel has been dodified to speed up the program
3862// - Some minor changes have been done to remove warnings
3863// - buffer size and split version of the Branches have been removed
3864// - Some modifications were needed to adat the program to the new
3865// MRawEvtData::DeletePixels
3866//
3867// Revision 1.28 2001/10/26 16:31:45 blanch
3868// Removing several warnings.
3869//
3870// Revision 1.27 2001/09/05 10:04:33 blanch
3871// *** empty log message ***
3872//
3873// Revision 1.26 2001/07/19 09:29:53 blanch
3874// Different threshold for each pixel can be used.
3875//
3876// Revision 1.25 2001/05/08 08:07:54 blanch
3877// New numbering for branches from different trigger conditions has been
3878// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
3879// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
3880// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
3881// information already filled.
3882//
3883// Revision 1.24 2001/04/06 16:48:09 magicsol
3884// New camera version able to read the new format of the reflector output:
3885// reflector 0.4
3886//
3887// Revision 1.23 2001/03/28 16:13:41 blanch
3888// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
3889// been solved.
3890//
3891// Revision 1.22 2001/03/20 18:52:43 blanch
3892// *** empty log message ***
3893//
3894// Revision 1.21 2001/03/19 19:53:03 blanch
3895// Some print outs have been removed.
3896//
3897// Revision 1.20 2001/03/19 19:30:06 magicsol
3898// Minor changes have been done to improve the FADC pedestals treatment.
3899//
3900// Revision 1.19 2001/03/05 11:14:41 magicsol
3901// I changed the position of readinf a parameter. It is a minnor change.
3902//
3903// Revision 1.18 2001/03/05 10:36:52 blanch
3904// A branch with information about the FADC simulation (MMcFadcHeader) is writen
3905// in the McHeaders tree of the aoutput root file.
3906// The NSB contribution is only applied if the the number of phe form the shower
3907// are above a given number.
3908//
3909// Revision 1.17 2001/02/23 11:05:57 magicsol
3910// Small changes due to slightly different output format and the introduction of
3911// pedesals for teh FADC.
3912//
3913// Revision 1.16 2001/01/18 18:44:40 magicsol
3914// *** empty log message ***
3915//
3916// Revision 1.15 2001/01/17 09:32:27 harald
3917// The changes are neccessary to have the same name for trees in MC and in
3918// data. So now there should be no differences in MC data and real data from
3919// FADC system.
3920//
3921// Revision 1.14 2001/01/15 12:33:34 magicsol
3922// Some modifications have been done to use the new (Dec'2000) Raw data format.
3923// There are also some minnors modifications to adapt the improvements in the
3924// MTrigger class (overlaping time and trigger cells).
3925//
3926// Revision 1.13 2000/10/25 08:14:23 magicsol
3927// The routine that produce poisson random numbers to decide how many phe
3928// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
3929// is used.
3930//
3931// Revision 1.12 2000/09/22 17:40:18 harald
3932// Added a lot of changes done by oscar.
3933//
3934// Revision 1.11 2000/09/21 11:47:33 harald
3935// Oscar found some smaller errors in the calculation of the pixel shape and
3936// corrected it.
3937//
3938// $Log: not supported by cvs2svn $
3939// Revision 1.46 2002/11/08 17:51:00 blanch
3940// *** empty log message ***
3941//
3942// Revision 1.45 2002/10/29 17:15:27 blanch
3943// Reading several reflector versions.
3944//
3945// Revision 1.44 2002/10/18 16:53:03 blanch
3946// Modification to read several reflector files.
3947//
3948// Revision 1.43 2002/09/13 10:53:39 blanch
3949// Minor change to remove some undisired comments.
3950//
3951// Revision 1.42 2002/09/09 16:00:49 blanch
3952// Statement has been included to avoid writting to disk MParContainer and MArray.
3953// It has also been added the effect of the WC, the actual values must be added,
3954// once they are measured.
3955//
3956// Revision 1.41 2002/09/04 09:57:42 blanch
3957// Modifications done to use MGeomCam from MARS.
3958//
3959// Revision 1.40 2002/07/16 16:15:22 blanch
3960// A first implementation of the Star field rotation has been done.
3961//
3962// Revision 1.39 2002/04/27 10:48:39 blanch
3963// Some unused varibles have been removed.
3964//
3965// Revision 1.38 2002/03/18 18:44:29 blanch
3966// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
3967//
3968// Revision 1.37 2002/03/18 16:42:20 blanch
3969// The data member fProductionSite of the MMcRunHeader has been set to 0,
3970// which means that the prodution site is unknown.
3971//
3972// Revision 1.35 2002/03/15 15:15:52 blanch
3973// Several modifications to simulate the actual trigger zone.
3974//
3975// Revision 1.34 2002/03/13 18:13:56 blanch
3976// Some changes to fill correctly the new format of MMcRunHeader.
3977//
3978// Revision 1.33 2002/03/04 17:21:48 blanch
3979// Small and not important changes.
3980//
3981// Revision 1.32 2002/02/28 15:04:52 blanch
3982// A small back has been solved. Before, while not using the option
3983// writte_all_images, not all triggered showers were stored. Now it is solved.
3984// For that it is importatn taht the less restrictive trigger option is
3985// checked first.
3986// A new facility has been introduced and now one can choose the step size in
3987// trigger loop mode for the discriminator threshold.
3988// The close-compact topology for less than 3 pixels does not make sense. Before
3989// the program was ignoring that, now it switch to simple neighbour condition.
3990//
3991// Revision 1.31 2002/01/18 17:41:02 blanch
3992// The option of adding noise to all pixels or to not adding the noise
3993// has been added.
3994// We removed the pixels larger than 577. When there were more than one
3995// trigger in one shower, the pixel number was increasing. Now it is
3996// flagged by the variable MMcTrig::fFirstLvlTrig.
3997//
3998// Revision 1.30 2001/11/27 09:49:54 blanch
3999// Fixing bug which was treating wrongly the extension of star photons.
4000//
4001// Revision 1.29 2001/11/14 17:38:23 blanch
4002// Sveral changes have been done:
4003// - bpoint_in_in_pixel has been dodified to speed up the program
4004// - Some minor changes have been done to remove warnings
4005// - buffer size and split version of the Branches have been removed
4006// - Some modifications were needed to adat the program to the new
4007// MRawEvtData::DeletePixels
4008//
4009// Revision 1.28 2001/10/26 16:31:45 blanch
4010// Removing several warnings.
4011//
4012// Revision 1.27 2001/09/05 10:04:33 blanch
4013// *** empty log message ***
4014//
4015// Revision 1.26 2001/07/19 09:29:53 blanch
4016// Different threshold for each pixel can be used.
4017//
4018// Revision 1.25 2001/05/08 08:07:54 blanch
4019// New numbering for branches from different trigger conditions has been
4020// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
4021// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
4022// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
4023// information already filled.
4024//
4025// Revision 1.24 2001/04/06 16:48:09 magicsol
4026// New camera version able to read the new format of the reflector output:
4027// reflector 0.4
4028//
4029// Revision 1.23 2001/03/28 16:13:41 blanch
4030// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
4031// been solved.
4032//
4033// Revision 1.22 2001/03/20 18:52:43 blanch
4034// *** empty log message ***
4035//
4036// Revision 1.21 2001/03/19 19:53:03 blanch
4037// Some print outs have been removed.
4038//
4039// Revision 1.20 2001/03/19 19:30:06 magicsol
4040// Minor changes have been done to improve the FADC pedestals treatment.
4041//
4042// Revision 1.19 2001/03/05 11:14:41 magicsol
4043// I changed the position of readinf a parameter. It is a minnor change.
4044//
4045// Revision 1.18 2001/03/05 10:36:52 blanch
4046// A branch with information about the FADC simulation (MMcFadcHeader) is writen
4047// in the McHeaders tree of the aoutput root file.
4048// The NSB contribution is only applied if the the number of phe form the shower
4049// are above a given number.
4050//
4051// Revision 1.17 2001/02/23 11:05:57 magicsol
4052// Small changes due to slightly different output format and the introduction of
4053// pedesals for teh FADC.
4054//
4055// Revision 1.16 2001/01/18 18:44:40 magicsol
4056// *** empty log message ***
4057//
4058// Revision 1.15 2001/01/17 09:32:27 harald
4059// The changes are neccessary to have the same name for trees in MC and in
4060// data. So now there should be no differences in MC data and real data from
4061// FADC system.
4062//
4063// Revision 1.14 2001/01/15 12:33:34 magicsol
4064// Some modifications have been done to use the new (Dec'2000) Raw data format.
4065// There are also some minnors modifications to adapt the improvements in the
4066// MTrigger class (overlaping time and trigger cells).
4067//
4068// Revision 1.13 2000/10/25 08:14:23 magicsol
4069// The routine that produce poisson random numbers to decide how many phe
4070// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
4071// is used.
4072//
4073// Revision 1.10 2000/07/04 14:10:20 MagicSol
4074// Some changes have been done in the root output file. The RawEvt tree is only
4075// stored in the single trigger mode.
4076// The trigger input parameters are also given by the general input card.
4077// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
4078// each one can be studied seperately.
4079//
4080// Revision 1.9 2000/06/13 13:25:24 blanch
4081// The multiple arrays have been replaced, since they do not work
4082// in alpha machines. Now we are using pointers and the command new
4083// to allocate memory.
4084//
4085// Revision 1.8 2000/05/11 13:57:27 blanch
4086// The option to loop over several trigger configurations has been included.
4087// Some non-sense with triggertime range has been solved.
4088// Montecarlo information and ADC counts are saved in a root file.
4089// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
4090//
4091// Revision 1.7 2000/03/24 18:10:46 blanch
4092// A first FADC simulation and a trigger simulation are already implemented.
4093// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
4094// A loop over trigger threshold and some corretcions in the time range where it looks for a trigger will be implemented as soon as possible.
4095//
4096// Revision 1.6 2000/03/20 18:35:11 blanch
4097// The trigger is already implemented but it does not save the trigger information in any file as it is implemented in timecam. In the next days there will be a version which also creates the files with the trigger information. It is going to be a mixing of the current camera and timecam programs.
4098//
4099// Revision 1.5 2000/02/18 17:40:35 petry
4100// This version includes drastic changes compared to camera.cxx 1.4.
4101// It is not yet finished and not immediately useful because the
4102// trigger simulation is not yet re-implemented. I had to take it
4103// out together with some other stuff in order to tidy the whole
4104// program up. This is not meant as an insult to anyone. I needed
4105// to do this in order to be able to work on it.
4106//
4107// This version has been put in the repository in order to be
4108// able to share the further development with others.
4109//
4110// If you need something working, wait or take an earlier one.
4111// See file README.
4112//
4113// Revision 1.4 2000/01/25 08:36:23 petry
4114// The pixelization in previous versions was buggy.
4115// This is the first version with a correct pixelization.
4116//
4117// Revision 1.3 2000/01/20 18:22:17 petry
4118// Found little bug which makes camera crash if it finds a photon
4119// of invalid wavelength. This bug is now fixed and the range
4120// of valid wavelengths extended to 290 - 800 nm.
4121// This is in preparation for the NSB simulation to come.
4122// Dirk
4123//
4124// Revision 1.2 1999/11/19 08:40:42 harald
4125// Now it is possible to compile the camera programm under osf1.
4126//
4127// Revision 1.1.1.1 1999/11/05 11:59:31 harald
4128// This the starting point for CVS controlled further developments of the
4129// camera program. The program was originally written by Jose Carlos.
4130// But here you can find a "rootified" version to the program. This means
4131// that there is no hbook stuff in it now. Also the output of the
4132// program changed to the MagicRawDataFormat.
4133//
4134// The "rootification" was done by Dirk Petry and Harald Kornmayer.
4135//
4136// Revision 1.3 1999/10/22 15:01:28 petry
4137// version sent to H.K. and N.M. on Fri Oct 22 1999
4138//
4139// Revision 1.2 1999/10/22 09:44:23 petry
4140// first synthesized version which compiles and runs without crashing;
4141//
4142// Revision 1.1.1.1 1999/10/21 16:35:10 petry
4143// first synthesised version
4144//
4145// Revision 1.13 1999/03/15 14:59:05 gonzalez
4146// camera-1_1
4147//
4148// Revision 1.12 1999/03/02 09:56:10 gonzalez
4149// *** empty log message ***
4150//
4151//
4152//!@}
4153
4154//=EOF
Note: See TracBrowser for help on using the repository browser.