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

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