source: trunk/MagicSoft/Simulation/Detector/TimeCam/timecam.cxx

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