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

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