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

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