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

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