source: branches/start/MagicSoft/Simulation/Detector/TimeCam/timecam.cxx@ 20115

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