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

Last change on this file since 1089 was 1089, checked in by blanch, 23 years ago
Fixing bug which was treating wrongly the extension of star photons.
File size: 98.1 KB
Line 
1//!/////////////////////////////////////////////////////////////////////
2//
3// camera
4//
5// @file camera.cxx
6// @title Camera simulation
7// @subtitle Code for the simulation of the camera phase
8// @desc Code for the simulation of the camera of CT1 and MAGIC
9// @author J C Gonzalez
10// @email gonzalez@mppmu.mpg.de
11// @date Thu May 7 16:24:22 1998
12//
13//----------------------------------------------------------------------
14//
15// Created: Thu May 7 16:24:22 1998
16// Author: Jose Carlos Gonzalez
17// Purpose: Program for reflector simulation
18// Notes: See files README for details
19//
20//----------------------------------------------------------------------
21//
22// $RCSfile: camera.cxx,v $
23// $Revision: 1.30 $
24// $Author: blanch $
25// $Date: 2001-11-27 09:49:54 $
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 "TRandom.h"
52#include "TApplication.h"
53
54#include "TFile.h"
55#include "TTree.h"
56#include "TBranch.h"
57#include "TCanvas.h"
58
59#include "TArrayC.h"
60
61#include "MTrigger.hxx"
62#include "MFadc.hxx"
63#include "MLons.hxx"
64
65#include "MRawRunHeader.h"
66#include "MRawEvtData.h"
67#include "MRawEvtHeader.h"
68#include "MRawCrateData.h"
69#include "MMcEvt.hxx"
70#include "MMcTrig.hxx"
71#include "MMcRunHeader.hxx"
72#include "MMcTrigHeader.hxx"
73#include "MMcFadcHeader.hxx"
74
75/*!@"
76
77 All the defines are located in the file |camera.h|.
78
79 @"*/
80
81#include "camera.h"
82
83//!@}
84
85/*!@"
86
87 The following set of flags are used in time of compilation. They do
88 not affect directly the behaviour of the program at run-time
89 (though, of course, if you disconnected the option for
90 implementation of the Trigger logic, you will not be able to use any
91 trigger at all. The 'default' values mean default in the sense of
92 what you got from the server when you obtained this program.
93
94 @"*/
95
96//!@{
97
98// flag for debugging (default: OFF )
99#define __DEBUG__
100#undef __DEBUG__
101
102
103//!@}
104
105//=-----------------------------------------------------------
106//!@subsection Definition of global variables.
107
108/*!@"
109
110 Now we define some global variables with data about the telescope,
111 such as "focal distance", number of pixels/mirrors,
112 "size of the camera", and so on.
113
114 @"*/
115
116/*!@"
117
118 Depending on the telescope we are using (CT1 or MAGIC), the
119 information stored in the definition file is different.
120 The variable |ct_Type| has the value 0 when we use
121 CT1, and 1 when we use MAGIC.
122
123 @"*/
124
125//!@{
126static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC
127//!@}
128
129/*!@"
130
131 And this is the information about the whole telescope.
132
133 @"*/
134
135//!@{
136
137// parameters of the CT (from the CT definition file)
138
139////@: Focal distances [cm]
140//static float *ct_Focal;
141
142//@: Mean Focal distances [cm]
143static float ct_Focal_mean;
144
145//@: STDev. Focal distances [cm]
146static float ct_Focal_std;
147
148//@: Mean Point Spread function [cm]
149static float ct_PSpread_mean;
150
151//@: STDev. Point Spread function [cm]
152static float ct_PSpread_std;
153
154//@: STDev. Adjustmente deviation [cm]
155static float ct_Adjustment_std;
156
157//@: Radius of the Black Spot in mirror [cm]
158static float ct_BlackSpot_rad;
159
160//@: Radius of one mirror [cm]
161static float ct_RMirror;
162
163//@: Camera width [cm]
164static float ct_CameraWidth;
165
166//@: Pixel width [cm]
167static float ct_PixelWidth;
168
169//@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)
170static float ct_PixelWidth_corner_2_corner;
171
172//@: ct_PixelWidth_corner_2_corner / 2
173static float ct_PixelWidth_corner_2_corner_half;
174
175//@: Number of mirrors
176static int ct_NMirrors = 0;
177
178//@: Number of pixels
179static int ct_NPixels;
180
181//@: Number of pixels
182static int ct_NCentralPixels;
183
184//@: Number of pixels
185static int ct_NGapPixels;
186
187//@: ct_Apot = ct_PixelWidth / 2
188static float ct_Apot;
189
190//@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth
191static float ct_2Apot;
192
193//@: name of the CT definition file to use
194static char ct_filename[256];
195
196//@: list of showers to be skipped
197static int *Skip;
198
199//@: number of showers to be skipped
200static int nSkip=0;
201
202//@: flag: TRUE: data come from STDIN; FALSE: from file
203static int Data_From_STDIN = FALSE;
204
205//@: flag: TRUE: write all images to output; FALSE: only triggered showers
206static int Write_All_Images = FALSE;
207
208//@: flag: TRUE: write all data to output; FALSE: only triggered showers
209static int Write_All_Data = FALSE;
210
211static int Write_McEvt = TRUE;
212static int Write_McTrig = TRUE;
213static int Write_McFADC = TRUE;
214static int Write_RawEvt = FALSE;
215
216//@: flag: TRUE: selection on the energy
217static int Select_Energy = TRUE;
218
219//@: Lower edge of the selected energy range (in GeV)
220static float Select_Energy_le = 0.0;
221
222//@: Upper edge of the selected energy range (in GeV)
223static float Select_Energy_ue = 100000.0;
224
225//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
226static int FADC_Scan = FALSE;
227
228//@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't
229static int Trigger_Scan = FALSE;
230
231//@: flag: TRUE: loop trigger analysis over several thresholds, multiplicities and topologies; FALSE: a single trigger configuration
232static int Trigger_Loop = FALSE;
233
234//@: flag: TRUE: Different threshold for each pixel ; FALSE: same threshold for all pixels
235static int Individual_Thres_Pixel = FALSE;
236
237//@: Properties of the trigger
238static float Trigger_gate_length = 6.0;
239static float Trigger_response_ampl = 1.0;
240static float Trigger_response_fwhm = 2.0;
241static float Trigger_overlaping_time= 0.25;
242
243//@: Properties of the FADC
244static float FADC_response_ampl = MFADC_RESPONSE_AMPLITUDE;
245static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
246
247//@: Trigger conditions for a single trigger mode
248static float qThreshold[TRIGGER_PIXELS];
249static int Trigger_multiplicity = 4;
250static int Trigger_topology = 2;
251
252//@: Upper and lower edges of the trigger loop
253static int Trigger_loop_lthres = 0;
254static int Trigger_loop_uthres = 10;
255static int Trigger_loop_lmult = 2;
256static int Trigger_loop_umult = 10;
257static int Trigger_loop_ltop = 0;
258static int Trigger_loop_utop = 2;
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
275//@: coordinates x,y for each pixel
276static float **pixary;
277
278//@: indexes of pixels neighbours of a given one
279static int **pixneig;
280
281//@: number of neighbours a pixel have
282static int *npixneig;
283
284//@: contents of the pixels (ph.e.)
285static float *fnpix;
286
287
288//!@}
289
290/*!@"
291
292 The following double-pointer is a 2-dimensional table with the
293 Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
294 of the wavelength @$\lambda@$. The routine |read_pixels()| will read
295 also this information from the file |qe.dat|.
296
297 @"*/
298
299//!@{
300// Pointer to a table with QE, number of datapoints, and wavelengths
301
302//@: table of QE
303static float ***QE;
304
305//@: number of datapoints for the QE curve
306static int pointsQE;
307
308//@: table of QE
309static float *QElambda;
310//!@}
311
312/*!@"
313
314 The following double-pointer is a 2-dimensional table with information
315 about each mirror in the dish. The routine |read_ct_file()| will read
316 this information from the CT definition file.
317
318 @"*/
319
320//!@{
321// Pointer to a table with the following info.:
322
323static float **ct_data;
324
325/*
326 * TYPE=0 (CT1)
327 * i s rho theta x y z thetan phin xn yn zn
328 *
329 * i : number of the mirror
330 * s : arc length [cm]
331 * rho : polar rho of the position of the center of the mirror [cm]
332 * theta : polar angle of the position of the center of the mirror [cm]
333 * x : x coordinate of the center of the mirror [cm]
334 * y : y coordinate of the center of the mirror [cm]
335 * z : z coordinate of the center of the mirror [cm]
336 * thetan : polar theta angle of the direction where the mirror points to
337 * phin : polar phi angle of the direction where the mirror points to
338 * xn : xn coordinate of the normal vector in the center (normalized)
339 * yn : yn coordinate of the normal vector in the center (normalized)
340 * zn : zn coordinate of the normal vector in the center (normalized)
341 *
342 * TYPE=1 (MAGIC)
343 * i f sx sy x y z thetan phin
344 *
345 * i : number of the mirror
346 * f : focal distance of that mirror
347 * sx : curvilinear coordinate of mirror's center in X[cm]
348 * sy : curvilinear coordinate of mirror's center in X[cm]
349 * x : x coordinate of the center of the mirror [cm]
350 * y : y coordinate of the center of the mirror [cm]
351 * z : z coordinate of the center of the mirror [cm]
352 * thetan : polar theta angle of the direction where the mirror points to
353 * phin : polar phi angle of the direction where the mirror points to
354 * xn : xn coordinate of the normal vector in the center (normalized)
355 * yn : yn coordinate of the normal vector in the center (normalized)
356 * zn : zn coordinate of the normal vector in the center (normalized)
357 */
358//!@}
359
360/*!@"
361
362 We define a table into where random numbers will be stored.
363 The routines used for random number generation are provided by
364 |RANLIB| (taken from NETLIB, |www.netlib.org|), and by
365 the routine |double drand48(void)| (prototype defined in
366 |stdlib.h|) through the macro |RandomNumber| defined in
367 |camera.h|.
368
369 @"*/
370
371
372//!@}
373
374extern char FileName[];
375
376
377//=-----------------------------------------------------------
378// @subsection Main program.
379
380//!@{
381
382//++++++++++++++++++++++++++++++++++++++++
383// MAIN PROGRAM
384//----------------------------------------
385
386int main(int argc, char **argv)
387{
388
389 //!@' @#### Definition of variables.
390 //@'
391
392 char inname[256]; //@< input file name
393 char starfieldname[256]; //@< starfield input file name
394
395 char datname[256]; //@< data (ASCII) output file name
396
397 char rootname[256] ; //@< ROOT file name
398 char rootname_loop[256] ; //@< ROOT file name
399
400 char parname[256]; //@< parameters file name
401
402 char nsbpathname[256]; //@< directory with the dataabse for th NSB
403
404 char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
405
406 FILE *inputfile; //@< stream for the input file
407 ofstream datafile; //@< stream for the data file
408
409 MCEventHeader mcevth; //@< Event Header class (MC)
410 MCCphoton cphoton; //@< Cherenkov Photon class (MC)
411
412 int inumphe; //@< number of photoelectrons in an event from showers
413 float inumphensb; //@< number of photoelectrons in an event from nsb
414
415 float arrtmin_ns; //@ arrival time of the first photoelectron
416 float arrtmax_ns; //@ arrival time of the last photoelectron
417
418 float thetaCT, phiCT; //@< parameters of a given shower
419 float thetashw, phishw; //@< parameters of a given shower
420 float coreD, coreX, coreY; //@< core position
421 float impactD; //@< impact parameter
422 float l1, m1, n1; //@< auxiliary variables
423 float l2, m2, n2; //@< auxiliary variables
424 float num, den; //@< auxiliary variables
425
426 int nshow=0; //@< partial number of shower in a given run
427 int ntshow=0; //@< total number of showers
428 int ncph=0; //@< partial number of photons in a given run
429 int ntcph=0; //@< total number of photons
430
431 int i, j, k; //@< simple counters
432
433 int simulateNSB; //@< Will we simulate NSB?
434 int nphe2NSB; //@< From how many phe will we simulate NSB?
435 float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel)
436 float diffnsb_phepns[iMAXNUMPIX]; //@< diffuse NSB values for each pixel
437 //@< derived from meanNSB
438 float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS]; //@< non-diffuse nsb
439 //@< photoelectron rates
440 float nsb_phepns[iMAXNUMPIX]; //@< NSB values for each pixel
441
442 Float_t nsb_trigresp[TRIGGER_TIME_SLICES]; //@< array to write the trigger
443 //@< response from the database
444 Float_t nsb_fadcresp[(Int_t) SLICES_MFADC]; //@< array to write the fadc
445 //@< response from the database
446 Byte_t trigger_map[((Int_t)(TRIGGER_PIXELS/8))+1];//@< Pixels on when the
447 //@< camera triggers
448 Float_t fadc_elecnoise[CAMERA_PIXELS]; //@< Electronic niose for each pixel
449 Float_t fadc_pedestals[CAMERA_PIXELS]; //@< array for fadc pedestals values
450
451 float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
452 EXTWAVEBAND1,
453 EXTWAVEBAND2,
454 EXTWAVEBAND3,
455 EXTWAVEBAND4,
456 EXTWAVEBAND5
457 };
458 float zenfactor=1.0; // correction factor calculated from the extinction
459
460 float qTailCut; //@< Tail Cut value
461 int nIslandsCut; //@< Islands Cut value
462 int countIslands=FALSE; //@< Will we count the islands?
463 int anaPixels;
464
465 float fCorrection; //@< Factor to apply to pixel values (def. 1.)
466
467 int ntrigger = 0; //@< number of triggers in the whole file
468 int btrigger = 0; //@< trigger flag
469 int ithrescount; //@< counter for loop over threshold trigger
470 int imulticount; //@< counter for loop over multiplicity trigger
471 int itopocount; //@< counter for loop over topology trigger
472 int icontrigger; //@< number of trigger conditions to be analised
473 UShort_t numPix; //@< number of sets of fadc written counts
474
475 float fpixelthres[TRIGGER_PIXELS]; //@< Threshold values
476
477 TArrayC *fadcValues; //@< the analog Fadc signal for pixels
478
479 float plateScale_cm2deg; //@< plate scale (deg/cm)
480 float degTriggerZone; //@< trigger area in the camera (radius, in deg.)
481
482 int still_in_loop = FALSE;
483
484 //@< Variables to fill the McRunHeader
485 Int_t sfRaH = 5;
486 Int_t sfRaM = 34;
487 Int_t sfRaS = 32;
488 Int_t sfDeD = 22;
489 Int_t sfDeM = 00;
490 Int_t sfDeS = 52;
491 Float_t shthetamax = 0.0;
492 Float_t shthetamin = 0.0;
493 Float_t shphimax = 0.0;
494 Float_t shphimin = 0.0;
495 Float_t telestheta = 0.0;
496 Float_t telesphi = 0.0;
497 Float_t sofftheta = 0.0;
498 Float_t soffphi = 0.0;
499 UInt_t corsika = 5200 ;
500
501
502 struct camera cam; // structure holding the camera definition
503
504
505 //!@' @#### Definition of variables for |getopt()|.
506 //@'
507
508 int ch, errflg = 0; //@< used by getopt
509
510 //!@' @#### Initialisation of some variables
511 //@'
512
513 for(i=0;i<iMAXNUMPIX;i++){
514 nsb_phepns[i]=0;
515 for(j=0;j<iNUMWAVEBANDS;j++)
516 nsbrate_phepns[i][j]=0.0; //@< Starlight rates initialised at 0.0
517 }
518 /*!@'
519
520 @#### Beginning of the program.
521
522 We start with the main program. First we (could) make some
523 presentation, and follows the reading of the parameters file (now
524 from the |stdin|), the reading of the CT parameters file, and the
525 creation of the output file, where the processed data will be
526 stored.
527
528 */
529
530 //++
531 // START
532 //--
533
534 // make unbuffered output
535
536 cout.setf ( ios::stdio );
537
538 // parse command line options (see reflector.h)
539
540 parname[0] = '\0';
541
542 optarg = NULL;
543 while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
544 switch (ch) {
545 case 'f':
546 strcpy(parname, optarg);
547 break;
548 case 'h':
549 usage();
550 break;
551 default :
552 errflg++;
553 }
554
555 // show help if error
556
557 if ( errflg>0 )
558 usage();
559
560 // make some sort of presentation
561
562 present();
563
564 // read parameters file
565
566 if ( strlen(parname) < 1 )
567 readparam(NULL);
568 else
569 readparam(parname);
570
571 // read data from file or from STDIN?
572
573 Data_From_STDIN = get_data_from_stdin();
574
575 // write all images, even those without trigger?
576
577 Write_All_Images = get_write_all_images();
578
579 // write all data (i.e., ph.e.s in pixels)
580
581 Write_All_Data = get_write_all_data();
582
583 Write_McEvt = get_write_McEvt() ;
584 Write_McTrig = get_write_McTrig() ;
585 Write_McFADC = get_write_McFadc() ;
586 Write_RawEvt = get_write_RawEvt() ;
587
588 FADC_Scan = get_FADC_Scan();
589 Trigger_Scan = get_Trigger_Scan();
590
591 Individual_Thres_Pixel = get_indi_thres_pixel();
592
593 get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm);
594
595 get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);
596
597 Trigger_Loop = get_Trigger_Loop(&Trigger_loop_lthres, &Trigger_loop_uthres, &Trigger_loop_lmult, &Trigger_loop_umult, &Trigger_loop_ltop, &Trigger_loop_utop);
598
599 icontrigger =(Trigger_loop_uthres-Trigger_loop_lthres+1)*
600 (Trigger_loop_umult-Trigger_loop_lmult+1)*
601 (Trigger_loop_utop-Trigger_loop_ltop+1);
602
603 if (!Trigger_Loop){
604 get_Trigger_Single (qThreshold, &Trigger_multiplicity, &Trigger_topology);
605 icontrigger=1;
606 }
607 else
608 get_threshold(qThreshold);
609
610 // get filenames
611
612 strcpy( inname, get_input_filename() );
613 strcpy( starfieldname, get_starfield_filename() );
614 strcpy( datname, get_data_filename() );
615 strcpy( rootname, get_root_filename() );
616 strcpy( rootname_loop, get_loop_filename() );
617 strcpy( ct_filename, get_ct_filename() );
618 strcpy( nsbpathname, get_nsb_directory() );
619
620 // get different parameters of the simulation
621
622 qTailCut = get_tail_cut();
623 simulateNSB = get_nsb( &meanNSB, &nphe2NSB );
624 countIslands = get_islands_cut( &nIslandsCut );
625
626 // get selections on the parameters
627
628 Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
629
630 // log filenames information
631
632 log(SIGNATURE,
633 "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n",
634 "Filenames",
635 "In", inname,
636 "Stars", starfieldname,
637 "NSB database", nsbpathname,
638 "CT", ct_filename,
639 "Data", datname,
640 "ROOT", rootname
641 );
642
643 // log Trigger information
644
645 if (Trigger_Loop) {
646 log(SIGNATURE,
647 "%s:\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s\n",
648 "Trigger Loop mode",
649 "Threshold",Trigger_loop_lthres,Trigger_loop_uthres,
650 "Multiplicity",Trigger_loop_lmult,Trigger_loop_umult,
651 "Topology",Trigger_loop_ltop,Trigger_loop_utop,
652 rootname_loop);
653 }
654 else if (Individual_Thres_Pixel == FALSE){
655 log(SIGNATURE,
656 "%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n",
657 "Single Trigger mode",
658 "Threshold",qThreshold[0],
659 "Multiplicity",Trigger_multiplicity,
660 "Topology",Trigger_topology);
661 }
662 else{
663 log(SIGNATURE,
664 "%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n",
665 "Single Trigger mode",
666 "Threshold","Different for each pixel",
667 "Multiplicity",Trigger_multiplicity,
668 "Topology",Trigger_topology);
669 }
670 // log flags information
671
672 log(SIGNATURE,
673 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
674 "Flags",
675 "Data_From_STDIN", ONoff(Data_From_STDIN),
676 "Write_All_Images", ONoff(Write_All_Images),
677 "Write_All_Data", ONoff(Write_All_Data));
678
679 // log flags information
680
681 log(SIGNATURE,
682 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
683 "Root output",
684 "Write_McEvt", ONoff(Write_McEvt),
685 "Write_McTrig", ONoff(Write_McTrig),
686 "Write_McFADC", ONoff(Write_McFADC),
687 "Write_RawEvt", ONoff(Write_RawEvt));
688
689 // log parameters information
690
691 log(SIGNATURE,
692 "%s:\n\t%20s: %f\n\t%20s: %f %s\n",
693 "Parameters",
694 "t0 (Tail-cut)", qTailCut,
695 "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB));
696
697 // log selections
698
699 log(SIGNATURE,
700 "%s:\n\t%20s: %s (%f:%f)\n",
701 "Selections:",
702 "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
703
704 // Definition and initialization of array to save trigger statistics
705
706 int ***ntriggerloop;
707
708 ntriggerloop= new int ** [(int) (Trigger_loop_uthres+1-Trigger_loop_lthres)];
709 for (ithrescount=0;ithrescount<=Trigger_loop_uthres-Trigger_loop_lthres;ithrescount++){
710 ntriggerloop[ithrescount]= new int * [Trigger_loop_umult-Trigger_loop_lmult+1];
711 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
712 ntriggerloop[ithrescount][imulticount]= new int [Trigger_loop_utop-Trigger_loop_ltop+1];
713 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
714 ntriggerloop[ithrescount][imulticount][itopocount]=0;
715 }
716 }
717 }
718
719 // set all random numbers seeds
720
721 setall( get_seeds(0), get_seeds(1) );
722
723 // get list of showers to evt. skip
724
725 nSkip = get_nskip_showers();
726
727 if (nSkip > 0) {
728 Skip = new int[ nSkip ];
729 get_skip_showers( Skip );
730
731 log(SIGNATURE, "There are some showers to skip:\n");
732 for (i=0; i<nSkip; ++i)
733 log(SIGNATURE, "\tshower # %d\n", Skip[i]);
734 }
735
736
737 // read parameters from the ct.def file
738
739 read_ct_file();
740
741 // read camera setup
742
743 read_pixels(&cam);
744
745 Int_t Lev0, Lev1, Lev2 ;
746
747 fadcValues = new TArrayC(FADC_SLICES);
748
749 // number of pixels for parameters
750
751 anaPixels = get_ana_pixels();
752 anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
753
754 // initialise ROOT
755
756 TROOT simple("simple", "MAGIC Telescope Monte Carlo");
757
758 // initialise instance of Trigger and FADC classes
759
760 MTrigger Trigger(Trigger_gate_length, Trigger_overlaping_time,
761 Trigger_response_ampl, Trigger_response_fwhm); //@< A instance of the Class MTrigger
762
763 // Initialise McTrig information class if we want to save trigger informtion
764
765 MMcTrig **McTrig = NULL;
766 MMcTrigHeader **HeaderTrig = NULL;
767 MMcFadcHeader **HeaderFadc = NULL;
768
769 if (Write_McTrig){
770
771 McTrig = new MMcTrig * [icontrigger];
772
773 for (i=0;i<icontrigger;i++) {
774 McTrig[i] = new MMcTrig();
775 }
776
777 HeaderTrig = new MMcTrigHeader * [icontrigger];
778
779 for (i=0;i<icontrigger;i++) {
780 HeaderTrig[i] = new MMcTrigHeader();
781 }
782 }
783
784 if (Write_McFADC){
785
786 HeaderFadc = new MMcFadcHeader * [icontrigger];
787
788 for (i=0;i<icontrigger;i++) {
789 HeaderFadc[i] = new MMcFadcHeader();
790 }
791 }
792
793 MFadc fadc(FADC_response_ampl,FADC_response_fwhm) ; //@< A instance of the Class MFadc
794
795 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
796 //
797 // Set the FADC pedestals for that run
798 // Some modifications
799 // mut be done to simulate a more realistic distribution of the pedestals.
800 // This simualtion is done int the SetPedestals methode inside the
801 // class MFadc
802 // Currentlly a input_pedestal array is declared with the pedestals.
803 // Thy can also be set randomly following a flat distribution.
804 //
805 /////////////////////////////////////////////////////////////////////
806
807 Float_t input_pedestals[CAMERA_PIXELS];
808
809 for(i=0;i<CAMERA_PIXELS;i++)
810 input_pedestals[i]=10.5;
811
812 fadc.SetPedestals(input_pedestals);
813
814 // Prepare the raw data output
815
816 // Header Tree
817
818 MRawRunHeader *RunHeader= new MRawRunHeader();
819 MMcRunHeader *McRunHeader = new MMcRunHeader();
820
821 // Header branch
822
823 MRawEvtHeader **EvtHeader = NULL;
824
825 if (Write_RawEvt) {
826 EvtHeader = new MRawEvtHeader * [icontrigger];
827
828 for (i=0;i<icontrigger;i++) {
829 EvtHeader[i] = new MRawEvtHeader();
830 }
831 }
832
833 // Data branch
834
835 MRawEvtData **EvtData = NULL; // Data branch
836
837 if (Write_RawEvt) {
838 EvtData = new MRawEvtData * [icontrigger];
839
840 for (i=0;i<icontrigger;i++) {
841 EvtData[i] = new MRawEvtData();
842 EvtData[i]->Init(RunHeader); // We need thr RunHeader to read
843 // number of pixels
844 }
845 }
846
847 MMcEvt *McEvt = new MMcEvt ();
848
849 //
850 // initalize a temporal ROOT file
851 //
852
853 TFile outfile_temp ( rootname , "RECREATE" );
854
855 // create a Tree for the Header Event
856 TTree HeaderTree("RunHeaders","Headers of Run");
857
858 // define branches of Header Tree
859
860 char help[4];
861
862 HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
863 &RunHeader);
864
865 HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
866 &McRunHeader);
867
868 if(!Trigger_Loop && Write_McTrig){
869
870 HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader",
871 &HeaderTrig[0]);
872 }
873 if (Trigger_Loop && Write_McTrig){
874 for(char branchname[20],i=0;i<icontrigger;i++){
875
876 sprintf(help,"%i",i+1);
877 strcpy (branchname, "MMcTrigHeader;");
878 strcat (branchname, & help[0]);
879 strcat (branchname, ".");
880 HeaderTree.Branch(branchname,"MMcTrigHeader",
881 &HeaderTrig[i]);
882 }
883 }
884
885 if(!Trigger_Loop && Write_McFADC){
886
887 HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader",
888 &HeaderFadc[0]);
889 }
890 if (Trigger_Loop && Write_McFADC){
891 for(char branchname[20],i=0;i<icontrigger;i++){
892
893 sprintf(help,"%i",i+1);
894 strcpy (branchname, "MMcFadcHeader;");
895 strcat (branchname, & help[0]);
896 strcat (branchname, ".");
897 HeaderTree.Branch(branchname,"MMcFadcHeader",
898 &HeaderFadc[i]);
899 }
900 }
901
902 // Fill branches for MRawRunHeader
903
904 RunHeader->SetMagicNumber(kMagicNumber);
905 RunHeader->SetFormatVersion(2);
906 RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
907 RunHeader->SetRunType(256);
908 RunHeader->SetRunNumber(0);
909 RunHeader->SetNumSamples(0,FADC_SLICES);
910 RunHeader->SetNumCrates(1);
911 RunHeader->SetNumPixInCrate(CAMERA_PIXELS);
912
913
914 // Fill branches for MMcTrigHeader
915
916 if(!Trigger_Loop && Write_McTrig){
917
918 HeaderTrig[0]->SetTopology((Short_t) Trigger_topology);
919 HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity);
920 HeaderTrig[0]->SetThreshold(qThreshold);
921 HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
922 HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
923 HeaderTrig[0]->SetOverlap(Trigger_overlaping_time);
924 HeaderTrig[0]->SetGate(Trigger_gate_length);
925 }
926 if(Trigger_Loop && Write_McTrig){
927
928 for (int iconcount=0,ithrescount=0;ithrescount<=Trigger_loop_uthres-Trigger_loop_lthres;ithrescount++){
929 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
930 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
931 HeaderTrig[iconcount]->SetTopology((Short_t) itopocount+Trigger_loop_ltop);
932 HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
933 for(int i=0;i<TRIGGER_PIXELS;i++){
934 fpixelthres[i]=
935 ((Float_t)(ithrescount+Trigger_loop_lthres)>=qThreshold[i])?
936 (Float_t)(ithrescount+Trigger_loop_lthres):qThreshold[i];
937 }
938 HeaderTrig[iconcount]->SetThreshold( fpixelthres);
939 HeaderTrig[iconcount]->SetAmplitud(Trigger_response_ampl);
940 HeaderTrig[iconcount]->SetFwhm(Trigger_response_fwhm);
941 HeaderTrig[iconcount]->SetOverlap(Trigger_overlaping_time);
942 HeaderTrig[iconcount]->SetGate(Trigger_gate_length);
943 iconcount++;
944 }
945 }
946 }
947 }
948
949 // Fill branches for MMcFadcHeader
950
951 for(i=0;i<CAMERA_PIXELS;i++){
952 fadc_elecnoise[i]=2.0;
953 }
954 fadc.GetPedestals(&fadc_pedestals[0]);
955
956 if(!Trigger_Loop && Write_McFADC){
957
958 HeaderFadc[0]->SetShape(0.0);
959 HeaderFadc[0]->SetAmplitud(FADC_response_ampl);
960 HeaderFadc[0]->SetFwhm(FADC_response_fwhm);
961 HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],anaPixels);
962 HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
963 }
964 if(Trigger_Loop && Write_McFADC){
965
966 for (int iconcount=0,ithrescount=0;ithrescount<=Trigger_loop_uthres-Trigger_loop_lthres;ithrescount++){
967 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
968 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
969 HeaderFadc[iconcount]->SetShape(0.0);
970 HeaderFadc[iconcount]->SetAmplitud(Trigger_response_ampl);
971 HeaderFadc[iconcount]->SetFwhm(Trigger_response_fwhm);
972 HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0],anaPixels);
973 HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
974 iconcount++;
975 }
976 }
977 }
978 }
979
980 // Fill the Header Tree with the current leaves of each branch
981 // HeaderTree.Fill() ;
982
983
984 // create a Tree for the Event data stream
985 TTree EvtTree("Events","Normal Triggered Events");
986
987 if (Write_McEvt){
988
989 EvtTree.Branch("MMcEvt","MMcEvt",
990 &McEvt);
991 }
992
993 if(!Trigger_Loop){
994
995 if (Write_RawEvt){
996 EvtTree.Branch("MRawEvtHeader","MRawEvtHeader",
997 &EvtHeader[0]);
998 EvtTree.Branch("MRawEvtData","MRawEvtData",
999 &EvtData[0]);
1000 }
1001 if (Write_McTrig){
1002 EvtTree.Branch("MMcTrig","MMcTrig",
1003 &McTrig[0]);
1004 }
1005 }
1006 else{
1007 if (Write_McTrig){
1008 for(char branchname[10],i=0;i<icontrigger;i++){
1009
1010 sprintf(help,"%i",i+1);
1011 strcpy (branchname, "MMcTrig;");
1012 strcat (branchname, & help[0]);
1013 strcat (branchname, ".");
1014 EvtTree.Branch(branchname,"MMcTrig",
1015 &McTrig[i]);
1016 }
1017 }
1018 }
1019
1020 if (Trigger_Loop && Write_RawEvt){
1021 for(char branchname[15],i=0;i<icontrigger;i++){
1022
1023 sprintf(help,"%i",i+1);
1024 strcpy (branchname, "MRawEvtHeader;");
1025 strcat (branchname, & help[0]);
1026 strcat (branchname, ".");
1027 EvtTree.Branch(branchname,"MRawEvtHeader",
1028 &EvtHeader[i]);
1029 }
1030 for(char branchname[15],i=0;i<icontrigger;i++){
1031
1032 sprintf(help,"%i",i+1);
1033 strcpy (branchname, "MRawEvtData;");
1034 strcat (branchname, & help[0]);
1035 strcat (branchname, ".");
1036 EvtTree.Branch(branchname,"MRawEvtData",
1037 &EvtData[i]);
1038 }
1039 }
1040
1041 TApplication theAppTrigger("App", &argc, argv);
1042
1043 if(FADC_Scan){
1044 if (gROOT->IsBatch()) {
1045 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1046 // return 1;
1047 }
1048 }
1049 if(FADC_Scan){
1050 //TApplication theAppFadc("App", &argc, argv);
1051
1052 if (gROOT->IsBatch()) {
1053 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1054 // return 1;
1055 }
1056 }
1057
1058 // set plate scale (deg/cm) and trigger area (deg)
1059
1060 plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;
1061
1062 if ( ! get_trigger_radius( &degTriggerZone ) )
1063 degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
1064
1065 if ( ! get_correction( &fCorrection ) )
1066 fCorrection = 1.0;
1067
1068 // prepare the NSB simulation
1069
1070 // Instance of the Mlons class
1071 MLons lons(Trigger_response_ampl, Trigger_response_fwhm,
1072 MFADC_RESPONSE_AMPLITUDE, MFADC_RESPONSE_FWHM);
1073
1074 lons.SetPath(nsbpathname);
1075
1076 if( simulateNSB ){
1077
1078 //
1079 // Calculate the non-diffuse NSB photoelectron rates
1080 //
1081
1082 k = produce_nsbrates( starfieldname,
1083 &cam,
1084 nsbrate_phepns);
1085
1086 if (k != 0){
1087 cout << "Error when reading starfield... \nExiting.\n";
1088 exit(1);
1089 }
1090
1091 // calculate diffuse rate correcting for the pixel size
1092
1093 for(i=0; i<cam.inumpixels; i++){
1094 diffnsb_phepns[i] = meanNSB *
1095 cam.dpixsizefactor[i] * cam.dpixsizefactor[i];
1096 }
1097
1098 // calculate nsb rate including diffuse and starlight
1099 // we also include the extinction effect
1100 for(j=0;j<iNUMWAVEBANDS;j++){
1101 // calculate the effect of the atmospheric extinction
1102
1103 zenfactor = pow(10., -0.4 * ext[j] );
1104
1105 for(i=0; i<cam.inumpixels;i++){
1106 nsb_phepns[i]+=diffnsb_phepns[i]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[i][j];
1107 }
1108 }
1109 }
1110
1111 //
1112 // Read the reflector file with the Cherenkov data
1113 //
1114
1115 // select input file
1116
1117 if ( Data_From_STDIN ) {
1118
1119 inputfile = stdin;
1120
1121 }
1122 else{
1123
1124 log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
1125 inputfile = fopen( inname, "r" );
1126 if ( inputfile == NULL )
1127 error( SIGNATURE, "Cannot open input file: %s\n", inname );
1128
1129 }
1130
1131 // get signature, and check it
1132
1133 if(check_reflector_file( inputfile )==FALSE){
1134 exit(1);
1135 }
1136
1137 // open data file
1138
1139 log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
1140 datafile.open( datname );
1141
1142 if ( datafile.bad() )
1143 error( SIGNATURE, "Cannot open data file: %s\n", datname );
1144
1145 // initializes flag
1146
1147 strcpy( flag, " \0" );
1148
1149 // allocate space for PMTs numbers of pixels
1150
1151 fnpix = new float [ ct_NPixels ];
1152
1153
1154 //!@' @#### Main loop.
1155 //@'
1156
1157 // get flag
1158
1159 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1160
1161 // loop over the file
1162
1163 still_in_loop = TRUE;
1164
1165 while (
1166 ((! Data_From_STDIN) && ( !feof(inputfile) ))
1167 ||
1168 (Data_From_STDIN && still_in_loop)
1169 ) {
1170
1171 // reading .rfl files
1172 if(!isA( flag, FLAG_START_OF_RUN )){
1173
1174 // We exit
1175 //error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
1176 // We break the main loop
1177 cout<<"Warning: Expected start of run flag, but found:"<<flag<<endl;
1178 cout<<" We break the main loop"<<endl;
1179 break;
1180 }
1181 else { // found start of run
1182
1183 nshow=0;
1184
1185 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1186
1187 while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
1188 //
1189 // Clear Trigger and Fadc
1190 //
1191 Trigger.Reset() ;
1192 Trigger.ClearFirst();
1193 Trigger.ClearZero();
1194 fadc.Reset() ;
1195
1196 ++nshow;
1197 log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
1198
1199 // get MCEventHeader
1200 fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
1201
1202 // calculate core distance and impact parameter
1203
1204 coreD = mcevth.get_core(&coreX, &coreY);
1205
1206 // calculate impact parameter (shortest distance betwee the original
1207 // trajectory of the primary (assumed shower-axis) and the
1208 // direction where the telescope points to
1209 //
1210 // we use the following equation, given that the shower core position
1211 // is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
1212 // and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
1213 // and (l2,m2,n2)
1214 //
1215 // | |
1216 // | x1-x2 y1-y2 z1-z2 |
1217 // | |
1218 // + | l1 m1 n1 |
1219 // - | |
1220 // | l2 m2 n2 |
1221 // | |
1222 // dist = ------------------------------------ ( > 0 )
1223 // [ |l1 m1|2 |m1 n1|2 |n1 l1|2 ]1/2
1224 // [ | | + | | + | | ]
1225 // [ |l2 m2| |m2 n2| |n2 l2| ]
1226 //
1227 // playing a little bit, we get this reduced for in our case:
1228 //
1229 //
1230 // dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
1231 // [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
1232 // 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
1233
1234 // read the direction of the incoming shower
1235
1236 thetashw = mcevth.get_theta();
1237 phishw = mcevth.get_phi();
1238
1239 // calculate vector for shower
1240
1241 l1 = sin(thetashw)*cos(phishw);
1242 m1 = sin(thetashw)*sin(phishw);
1243 n1 = cos(thetashw);
1244
1245 // read the deviation of the telescope with respect to the shower
1246
1247 mcevth.get_deviations ( &thetaCT, &phiCT );
1248
1249 if ( (thetaCT == 0.) && (phiCT == 0.) ) {
1250
1251 // CT was looking to the source (both lines are parallel)
1252 // therefore, we calculate the impact parameter as the distance
1253 // between the CT axis and the core position
1254
1255 impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1256
1257 } else {
1258
1259 // the shower comes off-axis
1260
1261 // obtain with this the final direction of the CT
1262
1263 thetaCT += thetashw;
1264 phiCT += phishw;
1265
1266 // calculate vector for telescope
1267
1268 l2 = sin(thetaCT)*cos(phiCT);
1269 m2 = sin(thetaCT)*sin(phiCT);
1270 n2 = cos(thetaCT);
1271
1272 num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
1273 den = (SQR(l1*m2 - l2*m1) +
1274 SQR(m1*n2 - m2*n1) +
1275 SQR(n1*l2 - n2*l1));
1276 den = sqrt(den);
1277
1278 impactD = fabs(num)/den;
1279
1280 }
1281
1282 // energy cut
1283
1284 if ( Select_Energy ) {
1285 if (( mcevth.get_energy() < Select_Energy_le ) ||
1286 ( mcevth.get_energy() > Select_Energy_ue )) {
1287 log(SIGNATURE, "select_energy: shower rejected.\n");
1288 continue;
1289 }
1290 }
1291
1292 // Read first and last time and put inumphe to 0
1293
1294 mcevth.get_times(&arrtmin_ns,&arrtmax_ns);
1295 inumphe=0;
1296
1297 // read the photons and produce the photoelectrons
1298
1299 k = produce_phes( inputfile,
1300 &cam,
1301 WAVEBANDBOUND1,
1302 WAVEBANDBOUND6,
1303 &Trigger, // will be changed by the function!
1304 &fadc, // will be changed by the function!
1305 &inumphe, // important for later: the size of photoe[]
1306 fnpix, // will be changed by the function!
1307 &ncph, // will be changed by the function!
1308 &arrtmin_ns, // will be changed by the function!
1309 &arrtmax_ns // will be changed by the function!
1310 );
1311
1312 if( k != 0 ){ // non-zero returnvalue means error
1313 cout << "Exiting.\n";
1314 exit(1);
1315 }
1316
1317 // NSB simulation
1318
1319 if(simulateNSB && nphe2NSB<=inumphe){
1320
1321 // Fill trigger and fadc response in the trigger class from the database
1322 for(i=0;i<cam.inumpixels;i++)
1323 if(nsb_phepns[i]>0.0){
1324 k=lons.GetResponse(nsb_phepns[i],0.1,
1325 & nsb_trigresp[0],& nsb_fadcresp[0]);
1326 if(k==0){
1327 cout << "Exiting.\n";
1328 exit(1);
1329 }
1330 Trigger.AddNSB(i,nsb_trigresp);
1331 fadc.AddSignal(i,nsb_fadcresp);
1332 }
1333 }// end if(simulateNSB && nphe2NSB<=inumphe) ...
1334
1335 log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1336 ncph, ntcph);
1337
1338 ntcph += ncph;
1339
1340 // skip it ?
1341
1342 for ( i=0; i<nSkip; ++i ) {
1343 if (Skip[i] == (nshow+ntshow)) {
1344 i = -1;
1345 break;
1346 }
1347 }
1348
1349 // if after the previous loop, the exit value of i is -1
1350 // then the shower number is in the list of showers to be
1351 // skipped
1352
1353 if (i == -1) {
1354 log(SIGNATURE, "\t\tskipped!\n");
1355 continue;
1356 }
1357
1358 cout << "Total number of phes: " << inumphe<<" + ";
1359 inumphensb=0;
1360 for (i=0;i<cam.inumpixels;i++){
1361 inumphensb+=nsb_phepns[i]*TOTAL_TRIGGER_TIME;
1362 }
1363 cout<<inumphensb<<endl;
1364
1365 //++++++++++++++++++++++++++++++++++++++++++++++++++
1366 // at this point we have a camera full of
1367 // ph.e.s
1368 // we should first apply the trigger condition,
1369 // and if there's trigger, then clean the image,
1370 // calculate the islands statistics and the
1371 // other parameters of the image (Hillas' parameters
1372 // and so on).
1373 //--------------------------------------------------
1374
1375 // TRIGGER HERE
1376
1377 //
1378 // now the noise of the electronic
1379 // (preamps, optical transmission,..) is introduced.
1380 // This is done inside the class MTrigger by the method ElecNoise.
1381 //
1382 Trigger.ElecNoise() ;
1383
1384 fadc.ElecNoise() ;
1385
1386 // We should simulate the AC coupling behaviour:
1387 // For the FADC it is only done for the NSB while producinr
1388 // the StarResponse database.
1389 // For the trigger is done in the Trigger.Diskriminate(), which
1390 // is called later (it should be seperated to speed up the program
1391
1392 // now a shift in the fadc signal due to the pedestlas is
1393 // intorduced
1394 // This is done inside the class MFadc by the method Pedestals
1395 fadc.Pedestals();
1396
1397 // We study several trigger conditons
1398 if(Trigger_Loop){
1399 // Loop over trigger threshold
1400 for (int iconcount=0,ithrescount=Trigger_loop_lthres;ithrescount<=Trigger_loop_uthres;ithrescount++){
1401 for (i=0;i<TRIGGER_PIXELS;i++)
1402 fpixelthres[i]=
1403 ((Float_t)(ithrescount)>=qThreshold[i])?
1404 (Float_t)(ithrescount):qThreshold[i];
1405 Trigger.SetThreshold(fpixelthres);
1406
1407 Trigger.Diskriminate();
1408 //
1409 // look if in all the signals in the trigger signal branch
1410 // is a possible Trigger. Therefore we habe to diskriminate all
1411 // the simulated analog signals (Method Diskriminate in class
1412 // MTrigger). We look simultanously for the moments at which
1413 // there are more than TRIGGER_MULTI pixels above the
1414 // CHANNEL_THRESHOLD.
1415 //
1416
1417 // Set trigger flags to zero
1418 Lev1=Lev2=0;
1419 btrigger=0;
1420
1421 // loop over multiplicity of trigger configuration
1422 for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1423 Trigger.SetMultiplicity(imulticount);
1424 Trigger.ClearZero();
1425
1426 Lev0=(Short_t) Trigger.ZeroLevel();
1427 if (Lev0>0 || Write_All_Images || btrigger){
1428
1429 // loop over topologies
1430 for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1431 Lev1=Lev2=0;
1432
1433 if(itopocount==0 && imulticount>7) continue;
1434 if(itopocount==2 && imulticount<3) continue;
1435 Trigger.SetTopology(itopocount);
1436 Trigger.ClearFirst();
1437
1438 //
1439 // Start the First Level Trigger simulation
1440 //
1441 if(Lev0!=0)
1442 Lev1=Trigger.FirstLevel();
1443 if (Write_McTrig)
1444 McTrig[iconcount]->SetFirstLevel (Lev1);
1445 if(Lev1>0) {
1446 btrigger= 1;
1447 ntriggerloop[ithrescount-Trigger_loop_lthres][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++;
1448 }
1449
1450 if(Lev1==0 && (Write_All_Images || btrigger)){
1451 btrigger= 1;
1452 Lev1=1;
1453 }
1454 numPix=0;
1455 for (Int_t ii=0;ii<Lev1;ii++){
1456 if (Write_McTrig)
1457 McTrig[iconcount]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
1458 if (Write_McTrig){
1459 Trigger.GetMapDiskriminator(trigger_map);
1460 McTrig[iconcount]->SetMapPixels(trigger_map,ii);
1461 }
1462 //
1463 // fill inside the class fadc the member output
1464 //
1465
1466 fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
1467
1468 if( Write_RawEvt ){
1469 //
1470 // Fill the header of this event
1471 //
1472
1473 EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
1474 // fill pixel information
1475 for(int i=0;i<ct_NPixels;i++){
1476 for (j=0;j<FADC_SLICES;j++){
1477 fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
1478 }
1479 EvtData[iconcount]->AddPixel(i+1000*ii,fadcValues,0);
1480 }
1481 }
1482 }
1483 //
1484 // Increase counter of analised trigger conditions
1485 //
1486 iconcount++;
1487 }
1488 }
1489 else{
1490 break;
1491 }
1492 }
1493 if (!btrigger) break;
1494 }
1495 if (btrigger){
1496
1497 //
1498 // fill the MMcEvt with all information
1499 //
1500
1501 if (Write_McEvt) {
1502 McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1503 mcevth.get_energy(),
1504 mcevth.get_theta(),
1505 mcevth.get_phi(),
1506 mcevth.get_core(),
1507 mcevth.get_coreX(),
1508 mcevth.get_coreY(),
1509 impactD,
1510 (UInt_t)mcevth.get_CORSIKA(),
1511 (UInt_t)mcevth.get_AtmAbs(),
1512 (UInt_t)mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot(),
1513 (UInt_t) ncph,
1514 (UInt_t) inumphe,
1515 (UInt_t) inumphensb+inumphe);
1516 }
1517 // Fill the Tree with the current leaves of each branch
1518 i=EvtTree.Fill() ;
1519
1520 // Clear the branches
1521 if(Write_McTrig){
1522 for(i=0;i<icontrigger;i++){
1523 McTrig[i]->Clear() ;
1524 }
1525 }
1526 if( Write_RawEvt ){
1527 for(i=0;i<icontrigger;i++){
1528 EvtHeader[i]->Clear() ;
1529 EvtData[i]->DeletePixels();
1530 }
1531 }
1532 if (Write_McEvt)
1533 McEvt->Clear() ;
1534 }
1535 }
1536 // We study a single trigger condition
1537 else {
1538
1539 // Setting trigger conditions
1540 Trigger.SetMultiplicity(Trigger_multiplicity);
1541 Trigger.SetTopology(Trigger_topology);
1542 for (int i=0;i<TRIGGER_PIXELS;i++)
1543 fpixelthres[i]=qThreshold[i];
1544 Trigger.SetThreshold(fpixelthres);
1545
1546 Trigger.Diskriminate() ;
1547
1548 //
1549 // look if in all the signals in the trigger signal branch
1550 // is a possible Trigger. Therefore we habe to diskriminate all
1551 // the simulated analog signals (Method Diskriminate in class
1552 // MTrigger). We look simultanously for the moments at which
1553 // there are more than TRIGGER_MULTI pixels above the
1554 // CHANNEL_THRESHOLD.
1555 //
1556
1557 Lev0 = (Short_t) Trigger.ZeroLevel() ;
1558
1559 Lev1 = Lev2 = 0 ;
1560
1561 //
1562 // Start the First Level Trigger simulation
1563 //
1564
1565 if ( Lev0 > 0 || Write_All_Images) {
1566 Lev1= Trigger.FirstLevel();
1567 if (Write_McTrig)
1568 McTrig[0]->SetFirstLevel (Lev1);
1569 }
1570 if (Lev1>0){
1571 ++ntrigger;
1572 }
1573 if (Lev1==0 && Write_All_Images){
1574 Lev1=1;
1575 }
1576
1577 numPix=0;
1578 for(Int_t ii=0;ii<Lev1;ii++){
1579 // Loop over different level one triggers
1580
1581 //
1582 // fill inside class fadc the member output
1583 //
1584 fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
1585
1586 if (Write_McTrig)
1587 McTrig[0]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
1588
1589 if (Write_McTrig){
1590 Trigger.GetMapDiskriminator(trigger_map);
1591 McTrig[0]->SetMapPixels(trigger_map,ii);
1592 }
1593
1594 // Fill Evt information
1595
1596 if (Write_RawEvt){
1597
1598 //
1599 // Fill the header of this event
1600 //
1601
1602 EvtHeader[0]->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ;
1603
1604 // fill pixel information
1605
1606 for(int i=0;i<ct_NPixels;i++){
1607 for (j=0;j<FADC_SLICES;j++){
1608 fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
1609 }
1610 EvtData[0]->AddPixel(i+ii*1000,fadcValues,0);
1611 }
1612 }
1613 //
1614 // fill the MMcEvt with all information
1615 //
1616
1617 if (Write_McEvt){
1618 McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1619 mcevth.get_energy(),
1620 mcevth.get_theta(),
1621 mcevth.get_phi(),
1622 mcevth.get_core(),
1623 mcevth.get_coreX(),
1624 mcevth.get_coreY(),
1625 impactD,
1626 (UInt_t)mcevth.get_CORSIKA(),
1627 (UInt_t)mcevth.get_AtmAbs(),
1628 (UInt_t)mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot(),
1629 (UInt_t) ncph,
1630 (UInt_t) inumphe,
1631 (UInt_t) inumphensb+inumphe);
1632 }
1633 // We don not count photons out of the camera.
1634
1635
1636 //
1637 // write it out to the file outfile
1638 //
1639
1640 EvtTree.Fill() ;
1641
1642 //
1643 // if a first level trigger occurred, then
1644 // 1. do some other stuff (not implemented)
1645 // 2. start the gui tool
1646
1647 if(FADC_Scan){
1648 if ( Lev0 > 0 ) {
1649 fadc.ShowSignal( McEvt, (Float_t) 60. ) ;
1650 }
1651 }
1652
1653 if(Trigger_Scan){
1654 if ( Lev0 > 0 ) {
1655 Trigger.ShowSignal(McEvt) ;
1656 }
1657 }
1658
1659 // clear all
1660 if (Write_RawEvt) EvtHeader[0]->Clear() ;
1661 if (Write_RawEvt) EvtData[0]->DeletePixels();
1662 if (Write_McEvt) McEvt->Clear() ;
1663 }
1664 if (Write_McTrig) McTrig[0]->Clear() ;
1665 }
1666
1667#ifdef __DEBUG__
1668 printf("\n");
1669
1670 for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1671
1672 for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1673
1674 if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1675
1676 if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1677
1678 printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1679 (int)pixels[ici][icj][PIXNUM],
1680 pixels[ici][icj][PIXX],
1681 pixels[ici][icj][PIXY],
1682 fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1683
1684 }
1685 }
1686 }
1687 }
1688
1689 for (int i=0; i<ct_NPixels; ++i) {
1690 printf("%d (%d): ", i, npixneig[i]);
1691 for (j=0; j<npixneig[i]; ++i)
1692 printf(" %d", pixneig[i][j]);
1693 printf("\n");
1694 }
1695
1696#endif // __DEBUG__
1697
1698
1699 // look for the next event
1700
1701 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1702
1703 } // end while there is a next event
1704
1705 if( !isA( flag, FLAG_END_OF_RUN )){
1706 error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1707 }
1708 else { // found end of run
1709 ntshow += nshow;
1710 log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1711
1712 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1713
1714 if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1715 log(SIGNATURE, "End of file . . .\n");
1716 still_in_loop = FALSE;
1717
1718 if ((! Data_From_STDIN) && ( !feof(inputfile) )){
1719
1720 // we have concatenated input files.
1721 // get signature of the next part and check it.
1722
1723 if(check_reflector_file( inputfile )==FALSE){
1724 exit(1);
1725 }
1726
1727 }
1728
1729 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1730
1731 } // end if found end of file
1732 } // end if found end of run
1733
1734 } // end if else found start of run
1735 } // end big while loop
1736
1737 //<@ Finally we should fill th McRunHeader
1738
1739
1740 get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
1741 get_teles_axis(&telestheta,&telesphi);
1742 get_source_off(&sofftheta,&soffphi);
1743 mcevth.get_theta_range(&shthetamin, &shthetamax);
1744 mcevth.get_phi_range(&shphimin,&shphimax);
1745 corsika=get_corsika_ver();
1746 if(!Trigger_Loop) icontrigger=0;
1747 McRunHeader->Fill(icontrigger,
1748 !Write_All_Images,
1749 Write_McEvt,
1750 Write_McTrig,
1751 Write_McFADC,
1752 CAMERA_PIXELS,
1753 (UInt_t)ntshow,
1754 (UInt_t)ntrigger,
1755 sfRaH,
1756 sfRaM,
1757 sfRaS,
1758 sfDeD,
1759 sfDeM,
1760 sfDeS,
1761 meanNSB,
1762 telestheta,
1763 telesphi,
1764 sofftheta,
1765 soffphi,
1766 shthetamax,
1767 shthetamin,
1768 shphimax,
1769 shphimin,
1770 corsika,
1771 (UInt_t)(REFL_VERSION*100),
1772 (UInt_t)(VERSION*100));
1773
1774 // Fill the Header Tree with the current leaves of each branch
1775 HeaderTree.Fill() ;
1776
1777 //++
1778 // put the Event to the root file
1779 //--
1780
1781 outfile_temp.Write() ;
1782 outfile_temp.Close() ;
1783
1784
1785 // close input file
1786
1787 log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
1788 ntshow, ntcph );
1789 datafile<<ntshow<<" event(s), with a total of "<<ntcph<<" C.photons"<<endl;
1790 if (Trigger_Loop){
1791 log( SIGNATURE, "Trigger Mode. \n");
1792 log( SIGNATURE, "Fraction of triggers: \n");
1793 datafile<<"Fraction of triggers: "<<endl;
1794 for (ithrescount=Trigger_loop_lthres;ithrescount<=Trigger_loop_uthres;ithrescount++){
1795 for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1796 for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1797 log( SIGNATURE, "Thres %d, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
1798 ithrescount,imulticount,itopocount,((float)ntriggerloop[ithrescount-Trigger_loop_lthres][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0), ntriggerloop[ithrescount-Trigger_loop_lthres][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop], ntshow);
1799 datafile<<"Thres "<<ithrescount<<", Multi "<<imulticount<<", Topo"<<itopocount<<": ";
1800 datafile<<((float)ntriggerloop[ithrescount-Trigger_loop_lthres][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0)<<"% ("<<ntriggerloop[ithrescount-Trigger_loop_lthres][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]<<" out of "<<ntshow<<")"<<endl;
1801 }
1802 }
1803 }
1804 }
1805 else{
1806 log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
1807 ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
1808 datafile<<"Fraction of triggers: "<<((float)ntrigger) / ((float)ntshow) * 100.0<<" ("<<ntrigger<<" out of "<<ntshow<<" )"<<endl;
1809 }
1810
1811 // close files
1812
1813 log( SIGNATURE, "Closing files\n" );
1814
1815 if( ! Data_From_STDIN ){
1816 fclose( inputfile );
1817 }
1818 datafile.close();
1819
1820 // program finished
1821
1822 log( SIGNATURE, "Done.\n");
1823
1824 return( 0 );
1825}
1826//!@}
1827
1828// @T \newpage
1829
1830//!@subsection Functions definition.
1831
1832//!-----------------------------------------------------------
1833// @name present
1834//
1835// @desc Make some presentation
1836//
1837// @date Sat Jun 27 05:58:56 MET DST 1998
1838//------------------------------------------------------------
1839// @function
1840
1841//!@{
1842void
1843present(void)
1844{
1845 cout << "##################################################\n"
1846 << SIGNATURE << '\n' << '\n'
1847 << "Processor of the reflector output\n"
1848 << "J C Gonzalez, Jun 1998\n"
1849 << "##################################################\n\n"
1850 << flush ;
1851}
1852//!@}
1853
1854
1855//!-----------------------------------------------------------
1856// @name usage
1857//
1858// @desc show help
1859//
1860// @date Tue Dec 15 16:23:30 MET 1998
1861//------------------------------------------------------------
1862// @function
1863
1864//!@{
1865void
1866usage(void)
1867{
1868 present();
1869 cout << "\nusage ::\n\n"
1870 << "\t camera "
1871 << " [ -@ paramfile ] "
1872 << " [ -h ] "
1873 << "\n\n or \n\n"
1874 << "\t camera < paramfile"
1875 << "\n\n";
1876 exit(0);
1877}
1878//!@}
1879
1880
1881//!-----------------------------------------------------------
1882// @name log
1883//
1884// @desc function to send log information
1885//
1886// @var funct Name of the caller function
1887// @var fmt Format to be used (message)
1888// @var ... Other information to be shown
1889//
1890// @date Sat Jun 27 05:58:56 MET DST 1998
1891//------------------------------------------------------------
1892// @function
1893
1894//!@{
1895void
1896log(const char *funct, char *fmt, ...)
1897{
1898 va_list args;
1899
1900 // Display the name of the function that called error
1901 printf("[%s]: ", funct);
1902
1903 // Display the remainder of the message
1904 va_start(args, fmt);
1905 vprintf(fmt, args);
1906 va_end(args);
1907}
1908//!@}
1909
1910
1911//!-----------------------------------------------------------
1912// @name error
1913//
1914// @desc function to send an error message, and abort the program
1915//
1916// @var funct Name of the caller function
1917// @var fmt Format to be used (message)
1918// @var ... Other information to be shown
1919//
1920// @date Sat Jun 27 05:58:56 MET DST 1998
1921//------------------------------------------------------------
1922// @function
1923
1924//!@{
1925void
1926error(const char *funct, char *fmt, ...)
1927{
1928 va_list args;
1929
1930 // Display the name of the function that called error
1931 fprintf(stdout, "ERROR in %s: ", funct);
1932
1933 // Display the remainder of the message
1934 va_start(args, fmt);
1935 vfprintf(stdout, fmt, args);
1936 va_end(args);
1937
1938 perror(funct);
1939
1940 exit(1);
1941}
1942//!@}
1943
1944
1945//!-----------------------------------------------------------
1946// @name isA
1947//
1948// @desc returns TRUE(FALSE), if the flag is(is not) the given
1949//
1950// @var s1 String to be searched
1951// @var flag Flag to compare with string s1
1952// @return TRUE: both strings match; FALSE: oth.
1953//
1954// @date Wed Jul 8 15:25:39 MET DST 1998
1955//------------------------------------------------------------
1956// @function
1957
1958//!@{
1959int
1960isA( char * s1, const char * flag ) {
1961 return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
1962}
1963//!@}
1964
1965
1966//!-----------------------------------------------------------
1967// @name read_ct_file
1968//
1969// @desc read CT definition file
1970//
1971// @date Sat Jun 27 05:58:56 MET DST 1998
1972//------------------------------------------------------------
1973// @function
1974
1975//!@{
1976void
1977read_ct_file(void)
1978{
1979 char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
1980 char token[ITEM_MAX_LENGTH]; //@< a single token
1981 int i, j; //@< dummy counters
1982
1983 log( "read_ct_file", "start.\n" );
1984
1985 ifstream ctin ( ct_filename );
1986
1987 if ( ctin.bad() )
1988 error( "read_ct_file",
1989 "Cannot open CT def. file: %s\n", ct_filename );
1990
1991 // loop till the "end" directive is reached
1992
1993 while (!ctin.eof()) {
1994
1995 // get line from stdin
1996
1997 ctin.getline(line, LINE_MAX_LENGTH);
1998
1999 // look for each item at the beginning of the line
2000
2001 for (i=0; i<=define_mirrors; i++)
2002 if (strstr(line, CT_ITEM_NAMES[i]) == line)
2003 break;
2004
2005 // if it is not a valid line, just ignore it
2006
2007 if (i == define_mirrors+1)
2008 continue;
2009
2010 // case block for each directive
2011
2012 switch ( i ) {
2013
2014 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2015
2016 // get focal distance
2017
2018 sscanf(line, "%s %d", token, &ct_Type);
2019
2020 log( "read_ct_file", "<Type of Telescope>: %s\n",
2021 ((ct_Type==0) ? "CT1" : "MAGIC") );
2022
2023 break;
2024
2025 case focal_distance: // <focal distance> [cm]
2026
2027 // get focal distance
2028
2029 sscanf(line, "%s %f", token, &ct_Focal_mean);
2030
2031 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2032
2033 break;
2034
2035 case focal_std: // s(focal distance) [cm]
2036
2037 // get focal distance
2038
2039 sscanf(line, "%s %f", token, &ct_Focal_std);
2040
2041 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2042
2043 break;
2044
2045 case point_spread: // <point spread> [cm]
2046
2047 // get point spread
2048
2049 sscanf(line, "%s %f", token, &ct_PSpread_mean);
2050
2051 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2052
2053 break;
2054
2055 case point_std: // s(point spread) [cm]
2056
2057 // get point spread
2058
2059 sscanf(line, "%s %f", token, &ct_PSpread_std);
2060
2061 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2062
2063 break;
2064
2065 case adjustment_dev: // s(adjustment_dev) [cm]
2066
2067 // get point spread
2068
2069 sscanf(line, "%s %f", token, &ct_Adjustment_std);
2070
2071 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2072
2073 break;
2074
2075 case black_spot: // radius of the black spot in the center [cm]
2076
2077 // get black spot radius
2078
2079 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2080
2081 log( "read_ct_file", "Radius of the black spots: %f cm\n",
2082 ct_BlackSpot_rad);
2083
2084 break;
2085
2086 case r_mirror: // radius of the mirrors [cm]
2087
2088 // get radius of mirror
2089
2090 sscanf(line, "%s %f", token, &ct_RMirror);
2091
2092 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2093
2094 break;
2095
2096 case n_mirrors: // number of mirrors
2097
2098 // get the name of the output_file from the line
2099
2100 sscanf(line, "%s %d", token, &ct_NMirrors);
2101
2102 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2103
2104 break;
2105
2106 case camera_width: // camera width [cm]
2107
2108 // get the name of the ct_file from the line
2109
2110 sscanf(line, "%s %f", token, &ct_CameraWidth);
2111
2112 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2113
2114 break;
2115
2116 case n_pixels: // number of pixels
2117
2118 // get the name of the output_file from the line
2119
2120 sscanf(line, "%s %d", token, &ct_NPixels);
2121
2122 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2123
2124 break;
2125
2126 case n_centralpixels: // number of central pixels
2127
2128 // get the name of the output_file from the line
2129
2130 sscanf(line, "%s %d", token, &ct_NCentralPixels);
2131
2132 log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
2133
2134 break;
2135
2136 case n_gappixels: // number of gap pixels
2137
2138 // get the name of the output_file from the line
2139
2140 sscanf(line, "%s %d", token, &ct_NGapPixels);
2141
2142 log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
2143
2144 break;
2145
2146 case pixel_width: // pixel width [cm]
2147
2148 // get the name of the ct_file from the line
2149
2150 sscanf(line, "%s %f", token, &ct_PixelWidth);
2151
2152 ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
2153 ct_PixelWidth_corner_2_corner_half =
2154 ct_PixelWidth_corner_2_corner * 0.50;
2155 ct_Apot = ct_PixelWidth / 2;
2156 ct_2Apot = ct_Apot * 2.0;
2157
2158 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2159
2160 break;
2161
2162 case define_mirrors: // read table with the parameters of the mirrors
2163
2164 log( "read_ct_file", "Table of mirrors data:\n" );
2165
2166 // check whether the number of mirrors was already set
2167
2168 if ( ct_NMirrors == 0 )
2169 error( "read_ct_file", "NMirrors was not set.\n" );
2170
2171 // allocate memory for paths list
2172
2173 log( "read_ct_file", "Allocating memory for ct_data\n" );
2174
2175 ct_data = new float*[ct_NMirrors];
2176
2177 for (i=0; i<ct_NMirrors; i++)
2178 ct_data[i] = new float[CT_NDATA];
2179
2180 // read data
2181
2182 log( "read_ct_file", "Reading mirrors data...\n" );
2183
2184 for (i=0; i<ct_NMirrors; i++)
2185 for (j=0; j<CT_NDATA; j++)
2186 ctin >> ct_data[i][j];
2187
2188 break;
2189
2190 } // switch ( i )
2191
2192 } // end while
2193
2194 // end
2195
2196 log( "read_ct_file", "done.\n" );
2197
2198 return;
2199}
2200//!@}
2201
2202
2203//!-----------------------------------------------------------
2204// @name read_pixels
2205//
2206// @desc read pixels data
2207//
2208// @date Fri Mar 12 16:33:34 MET 1999
2209//------------------------------------------------------------
2210// @function
2211
2212//!@{
2213void
2214read_pixels(struct camera *pcam)
2215{
2216 ifstream qefile;
2217 char line[LINE_MAX_LENGTH];
2218 int n, i, j, icount;
2219 float qe;
2220
2221 //------------------------------------------------------------
2222 // first, pixels' coordinates
2223
2224 pcam->inumpixels = ct_NPixels;
2225 pcam->inumcentralpixels = ct_NCentralPixels;
2226 pcam->inumgappixels = ct_NGapPixels;
2227 pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
2228 pcam->dpixdiameter_cm = ct_PixelWidth;
2229
2230 // initialize pixel numbers
2231
2232 pixary = new float* [2*ct_NCentralPixels];
2233 for ( i=0; i<2*ct_NCentralPixels; ++i )
2234 pixary[i] = new float[2];
2235
2236 pixneig = new int* [ct_NCentralPixels];
2237 for ( i=0; i<ct_NCentralPixels; ++i ) {
2238 pixneig[i] = new int[6];
2239 for ( j=0; j<6; ++j )
2240 pixneig[i][j] = -1;
2241 }
2242
2243 npixneig = new int[ct_NCentralPixels];
2244 for ( i=0; i<ct_NCentralPixels; ++i )
2245 npixneig[i] = 0;
2246
2247 // generate all coordinates
2248
2249 igen_pixel_coordinates(pcam);
2250
2251
2252 // calculate tables of neighbours
2253
2254#ifdef __DEBUG__
2255 for ( n=0 ; n<ct_NPixels ; ++n ) {
2256 cout << "Para el pixel " << n << ": ";
2257 for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2258 if ( pixels_are_neig(n,i) == TRUE ) {
2259 pixneig[n][npixneig[n]] = i;
2260 pixneig[i][npixneig[i]] = n;
2261 cout << i << ' ';
2262 ++npixneig[n];
2263 ++npixneig[i];
2264 }
2265 }
2266 cout << endl << flush;
2267 }
2268#else // ! __DEBUG__
2269 for ( n=0 ; n<ct_NCentralPixels ; ++n )
2270 for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2271 if ( pixels_are_neig(n,i) == TRUE ) {
2272 pixneig[n][npixneig[n]] = i;
2273 pixneig[i][npixneig[i]] = n;
2274 ++npixneig[n];
2275 ++npixneig[i];
2276 }
2277#endif // ! __DEBUG__
2278
2279#ifdef __DEBUG__
2280 for ( n=0 ; n<ct_NPixels ; ++n ) {
2281 cout << n << ':';
2282 for ( j=0; j<npixneig[n]; ++j)
2283 cout << ' ' << pixneig[n][j];
2284 cout << endl << flush;
2285 }
2286#endif // __DEBUG__
2287
2288 //------------------------------------------------------------
2289 // second, pixels' QE
2290
2291 // try to open the file
2292
2293 log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
2294
2295 qefile.open( QE_FILE );
2296
2297 // if it is wrong or does not exist, exit
2298
2299 if ( qefile.bad() )
2300 error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2301
2302 // read file
2303
2304 log("read_pixels", "Reading data . . .\n");
2305
2306 i=-1;
2307 icount = 0;
2308
2309 while ( ! qefile.eof() ) {
2310
2311 // get line from the file
2312
2313 qefile.getline(line, LINE_MAX_LENGTH);
2314
2315 // skip if comment
2316
2317 if ( *line == '#' )
2318 continue;
2319
2320 // if it is the first valid value, it is the number of QE data points
2321
2322 if ( i < 0 ) {
2323
2324 // get the number of datapoints
2325
2326 sscanf(line, "%d", &pointsQE);
2327
2328 // allocate memory for the table of QEs
2329
2330 QE = new float ** [ct_NPixels];
2331
2332 for ( i=0; i<ct_NPixels; ++i ) {
2333 QE[i] = new float * [2];
2334 QE[i][0] = new float[pointsQE];
2335 QE[i][1] = new float[pointsQE];
2336 }
2337
2338 QElambda = new float [pointsQE];
2339
2340 for ( i=0; i<pointsQE; ++i ) {
2341 qefile.getline(line, LINE_MAX_LENGTH);
2342 sscanf(line, "%f", &QElambda[i]);
2343 }
2344
2345 i=0;
2346
2347 continue;
2348 }
2349
2350 // get the values (num-pixel, num-datapoint, QE-value)
2351
2352 if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2353 break;
2354
2355 if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2356 ((j-1) < pointsQE) && ((j-1) > -1) ) {
2357 QE[i-1][0][j-1] = QElambda[j-1];
2358 QE[i-1][1][j-1] = qe;
2359 }
2360
2361 if ( i > ct_NPixels) break;
2362
2363 icount++;
2364
2365 }
2366
2367 if(icount/pointsQE < ct_NPixels){
2368 error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2369 icount/pointsQE, ct_NPixels );
2370 }
2371
2372 // close file
2373
2374 qefile.close();
2375
2376 // test QE
2377
2378 for(icount=0; icount< ct_NPixels; icount++){
2379 for(i=0; i<pointsQE; i++){
2380 if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
2381 QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
2382 error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2383 icount, i, QE[icount][0][i], QE[icount][1][i] );
2384 }
2385 }
2386 }
2387
2388 // end
2389
2390 log("read_pixels", "Done.\n");
2391
2392}
2393//!@}
2394
2395
2396//!-----------------------------------------------------------
2397// @name pixels_are_neig
2398//
2399// @desc check whether two pixels are neighbours
2400//
2401// @var pix1 Number of the first pixel
2402// @var pix2 Number of the second pixel
2403// @return TRUE: both pixels are neighbours; FALSE: oth.
2404//
2405// @date Wed Sep 9 17:58:37 MET DST 1998
2406//------------------------------------------------------------
2407// @function
2408
2409//!@{
2410int
2411pixels_are_neig(int pix1, int pix2)
2412{
2413 if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2414 SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2415 > ct_PixelWidth_corner_2_corner )
2416 return ( FALSE );
2417 else
2418 return ( TRUE );
2419}
2420//!@}
2421
2422//!-----------------------------------------------------------
2423// @name igen_pixel_coordinates
2424//
2425// @desc generate the pixel center coordinates
2426//
2427// @var *pcam structure camera containing all the
2428// camera information
2429// @return total number of pixels
2430//
2431// DP
2432//
2433// @date Thu Oct 14 10:41:03 CEST 1999
2434//------------------------------------------------------------
2435// @function
2436
2437//!@{
2438/******** igen_pixel_coordinates() *********************************/
2439
2440int igen_pixel_coordinates(struct camera *pcam) {
2441 /* generate pixel coordinates, return value is number of pixels */
2442
2443 int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2444 float fsegment_fract;
2445 double dtsize;
2446 double dhsize;
2447 double dpsize;
2448 double dxfirst_pix;
2449 double dyfirst_pix;
2450 double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
2451 double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
2452
2453
2454 double dstartx, dstarty; /* for the gap pixels and outer pixels */
2455 int j, nrow;
2456
2457 dpsize = pcam->dpixdiameter_cm;
2458 dtsize = dpsize * sqrt(3.) / 2.;
2459 dhsize = dpsize / 2.;
2460
2461 /* Loop over central pixels to generate co-ordinates */
2462
2463 for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2464
2465 /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2466
2467 pcam->dpixsizefactor[ipixno-1] = 1.;
2468
2469 in = 0;
2470
2471 i = 0;
2472 itot_inside_ring = 0;
2473 iring_no = 0;
2474
2475 /* Calculate the number of pixels out to and including the ring containing pixel number */
2476 /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2477
2478 while (itot_inside_ring == 0){
2479
2480 iN = 3*(i*(i+1)) + 1;
2481
2482 if (ipixno <= iN){
2483 iring_no = i;
2484 itot_inside_ring = iN;
2485 }
2486
2487 i++;
2488 }
2489
2490
2491 /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2492
2493 ipix_in_ring = 0;
2494 for (i = 0; i < iring_no; ++i){
2495
2496 ipix_in_ring = ipix_in_ring + 6;
2497 }
2498
2499 /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2500 /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2501 /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2502
2503 isegment = 0;
2504 fsegment_fract = 0.;
2505 if (iring_no > 0) {
2506
2507 isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2508
2509 fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2510
2511 }
2512
2513 /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2514 /* pixel width (flat to flat) dpsize. */
2515
2516 dxfirst_pix = dpsize*iring_no;
2517 dyfirst_pix = 0.;
2518
2519 /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2520
2521 ddxseg1 = - dhsize*iring_no;
2522 ddyseg1 = dtsize*iring_no;
2523 ddxseg2 = -dpsize*iring_no;
2524 ddyseg2 = 0.;
2525 ddxseg3 = ddxseg1;
2526 ddyseg3 = -ddyseg1;
2527 ddxseg4 = -ddxseg1;
2528 ddyseg4 = -ddyseg1;
2529 ddxseg5 = -ddxseg2;
2530 ddyseg5 = 0.;
2531 ddxseg6 = -ddxseg1;
2532 ddyseg6 = ddyseg1;
2533
2534 /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2535 /* anti-clockwise around the ring by adding the segment to segment vectors. */
2536
2537 switch (isegment) {
2538
2539 case 0:
2540
2541 pcam->dxc[ipixno-1] = 0.;
2542 pcam->dyc[ipixno-1] = 0.;
2543
2544 case 1:
2545 pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2546 pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2547
2548 break;
2549
2550 case 2:
2551
2552 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2553 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2554
2555 break;
2556
2557 case 3:
2558
2559 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
2560 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
2561
2562 break;
2563
2564 case 4:
2565
2566 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
2567 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
2568
2569 break;
2570
2571 case 5:
2572
2573 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
2574 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2575
2576 break;
2577
2578 case 6:
2579
2580 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2581 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2582
2583 break;
2584
2585 default:
2586
2587 fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2588 return(0);
2589
2590 } /* end switch */
2591
2592 } /* end for */
2593
2594 dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2595 dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2596
2597 if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2598
2599 j = pcam->inumcentralpixels;
2600
2601 for(i=0; i<pcam->inumgappixels; i=i+6){
2602 pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2603 pcam->dyc[j + i ] = dstarty;
2604 pcam->dpixsizefactor[j + i] = 1.;
2605 pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2606 pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2607 pcam->dpixsizefactor[j + i + 1] = 1.;
2608 pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2609 pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2610 pcam->dpixsizefactor[j + i+ 2] = 1.;
2611 pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2612 pcam->dyc[j + i + 3] = dstarty;
2613 pcam->dpixsizefactor[j + i+ 3] = 1.;
2614 pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2615 pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2616 pcam->dpixsizefactor[j + i+ 4] = 1.;
2617 pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2618 pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2619 pcam->dpixsizefactor[j + i + 5] = 1.;
2620 } /* end for */
2621 } /* end if */
2622
2623 /* generate positions of the outer pixels */
2624
2625 if( pcam->inumbigpixels > 0 ){
2626
2627 j = pcam->inumcentralpixels + pcam->inumgappixels;
2628
2629 for(i=0; i<pcam->inumbigpixels; i++){
2630 pcam->dpixsizefactor[j + i] = 2.;
2631 }
2632
2633 in = 0;
2634
2635 nrow = (int) ceil(dstartx / 2. / dpsize);
2636
2637 while(in < pcam->inumbigpixels){
2638
2639 pcam->dxc[j + in] = dstartx + dpsize;
2640 pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2641 pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2642 pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2643 pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2644 pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2645 pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2646 pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2647 pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2648 pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2649 pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2650 pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2651 for(i=1; i<nrow; i++){
2652 pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2653 pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2654 pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2655 pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2656 pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2657 pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2658 pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2659 pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2660 pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
2661 pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
2662 pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
2663 pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
2664 }
2665 in = in + 6 * nrow;
2666 dstartx = dstartx + 2. * dpsize;
2667 nrow = nrow + 1;
2668 } /* end while */
2669
2670 } /* end if */
2671
2672 return(pcam->inumpixels);
2673
2674}
2675//!@}
2676
2677//!-----------------------------------------------------------
2678// @name bpoint_is_in_pix
2679//
2680// @desc check if a point (x,y) in camera coordinates is inside a given pixel
2681//
2682// @var *pcam structure camera containing all the
2683// camera information
2684// @var dx, dy point coordinates in centimeters
2685// @var ipixnum pixel number (starting at 0)
2686// @return TRUE if the point is inside the pixel, FALSE otherwise
2687//
2688// DP
2689//
2690// @date Thu Oct 14 16:59:04 CEST 1999
2691//------------------------------------------------------------
2692// @function
2693
2694//!@{
2695
2696/******** bpoint_is_in_pix() ***************************************/
2697/*
2698int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
2699 // return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system)
2700 // the pixel is assumed to be a "closed set"
2701
2702 double a, b; // a = length of one of the edges of one pixel, b = half the width of one pixel
2703 double c, xx, yy; // auxiliary variable
2704
2705 b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
2706 a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
2707 c = 1./sqrt(3.);
2708
2709 if((ipixnum < 0)||(ipixnum >= pcam->inumpixels)){
2710 fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
2711 fprintf(stderr, "Exiting.\n");
2712 exit(203);
2713 }
2714 xx = dx - pcam->dxc[ipixnum];
2715 yy = dy - pcam->dyc[ipixnum];
2716
2717 if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) ||
2718 ((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
2719 return(TRUE); // point is inside
2720 }
2721 else{
2722 return(FALSE); // point is outside
2723 }
2724}
2725*/
2726
2727#define sqrt13 0.577350269 // = 1./sqrt(3.)
2728#define sqrt3 1.732050807 // = sqrt(3.)
2729
2730int bpoint_is_in_pix(double dx, double dy, struct camera *pcam)
2731{
2732 /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
2733 /* the pixel is assumed to be a "closed set" */
2734
2735 /*
2736 a = length of one of the edges of one pixel,
2737 b = half the width of one pixel
2738 */
2739
2740 const int num = pcam->inumpixels;
2741 const double diam = pcam->dpixdiameter_cm;
2742
2743 const double diam2 = diam/2;
2744 const double diam3 = diam/sqrt3;
2745
2746 for (int i=0; i<num; i++)
2747 {
2748 const double size = pcam->dpixsizefactor[i];
2749
2750 const double b = diam2 * size;
2751 const double a = diam3 * size;
2752
2753 const double xx = dx - pcam->dxc[i];
2754 const double yy = dy - pcam->dyc[i];
2755
2756 if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
2757 ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){
2758
2759 return i; // inside i
2760 }
2761
2762 // return -1; // outside
2763 }
2764 return -1; // outside
2765}
2766
2767//!@}
2768
2769//------------------------------------------------------------
2770// @name dist_r_P
2771//
2772// @desc distance straight line r - point P
2773//
2774// @date Sat Jun 27 05:58:56 MET DST 1998
2775// @function @code
2776//------------------------------------------------------------
2777// dist_r_P
2778//
2779// distance straight line r - point P
2780//------------------------------------------------------------
2781
2782float
2783dist_r_P(float a, float b, float c,
2784 float l, float m, float n,
2785 float x, float y, float z)
2786{
2787 return (
2788 sqrt((SQR((a-x)*m-(b-y)*l) +
2789 SQR((b-y)*n-(c-z)*m) +
2790 SQR((c-z)*l-(a-x)*n))/
2791 (SQR(l)+SQR(m)+SQR(n))
2792 )
2793 );
2794}
2795
2796//------------------------------------------------------------
2797// @name check_reflector_file
2798//
2799// @desc check if a given reflector file has the right signature
2800// @desc return TRUE or FALSE
2801//
2802// @date Mon Feb 14 16:44:21 CET 2000
2803// @function @code
2804//------------------------------------------------------------
2805
2806int check_reflector_file(FILE *infile){
2807
2808 char Signature[20]; // auxiliary variable
2809 char sign[20]; // auxiliary variable
2810
2811 strcpy(Signature, REFL_SIGNATURE);
2812
2813 strcpy(sign, Signature);
2814
2815 fread( (char *)sign, strlen(Signature), 1, infile);
2816
2817 if (strcmp(sign, Signature) != 0) {
2818 cout << "ERROR: Signature of .rfl file is not correct\n";
2819 cout << '"' << sign << '"' << '\n';
2820 cout << "should be: " << Signature << '\n';
2821 return(FALSE);
2822 }
2823
2824 fread( (char *)sign, 1, 1, infile);
2825
2826 return(TRUE);
2827
2828}
2829
2830//------------------------------------------------------------
2831// @name lin_interpol
2832//
2833// @desc interpolate linearly between two points returning the
2834// @desc y-value of the result
2835//
2836// @date Thu Feb 17 11:31:32 CET 2000
2837// @function @code
2838//------------------------------------------------------------
2839
2840float lin_interpol( float x1, float y1, float x2, float y2, float x){
2841
2842 if( (x2 - x1)<=0. ){ // avoid division by zero, return average
2843 cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
2844 return((y1+y2)/2.);
2845 }
2846 else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
2847 return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
2848 }
2849}
2850
2851
2852//------------------------------------------------------------
2853// @name produce_phes
2854//
2855// @desc read the photons of an event, pixelize them and simulate QE
2856// @desc return various statistics and the array of Photoelectrons
2857//
2858// @date Mon Feb 14 16:44:21 CET 2000
2859// @function @code
2860//------------------------------------------------------------
2861
2862int produce_phes( FILE *sp, // the input file
2863 struct camera *cam, // the camera layout
2864 float minwl_nm, // the minimum accepted wavelength
2865 float maxwl_nm, // the maximum accepted wavelength
2866 class MTrigger *trigger, // the generated phes
2867 class MFadc *fadc,
2868 int *itotnphe, // total number of produced photoelectrons
2869 float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
2870 int *incph, // total number of cph read
2871 float *tmin_ns, // minimum arrival time of all phes
2872 float *tmax_ns // maximum arrival time of all phes
2873 ){
2874
2875 static int i, k, ipixnum;
2876 static float cx, cy, wl, qe, t;
2877 static MCCphoton photon;
2878 static float **qept;
2879 static char flag[SIZE_OF_FLAGS + 1];
2880 static float radius;
2881
2882
2883 // reset variables
2884
2885 for ( i=0; i<cam->inumpixels; ++i ){
2886
2887 nphe[i] = 0.0;
2888
2889 }
2890
2891 *incph = 0;
2892
2893 radius = cam->dxc[cam->inumpixels-1]
2894 + 1.5*cam->dpixdiameter_cm*cam->dpixsizefactor[cam->inumpixels-1];
2895
2896 //- - - - - - - - - - - - - - - - - - - - - - - - -
2897 // read photons and "map" them into the pixels
2898 //--------------------------------------------------
2899
2900 // initialize CPhoton
2901
2902 photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
2903
2904 // read the photons data
2905
2906
2907 // loop over the photons
2908
2909 while (1) {
2910
2911 fread ( flag, SIZE_OF_FLAGS, 1, sp );
2912
2913 if (isA( flag, FLAG_END_OF_EVENT ))
2914 break;
2915
2916 memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
2917
2918 fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
2919
2920 // increase number of photons
2921
2922 (*incph)++;
2923
2924 // Chceck if photon is inside trigger time range
2925
2926 t = photon.get_t() ;
2927
2928 if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
2929 continue;
2930 }
2931
2932 //
2933 // Pixelization
2934 //
2935
2936 cx = photon.get_x();
2937 cy = photon.get_y();
2938
2939 // get wavelength
2940
2941 wl = photon.get_wl();
2942
2943 // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n";
2944
2945 // check if photon has valid wavelength and is inside outermost camera radius
2946
2947 if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){
2948 continue;
2949
2950 }
2951
2952 ipixnum = bpoint_is_in_pix(cx, cy, cam);
2953
2954 // -1 = the photon is in none of the pixels
2955 // 0 = the phton is in the central pixel, which is not used for trigger
2956 if (ipixnum==-1 || ipixnum==0) {
2957 continue;
2958 }
2959
2960 //+++
2961 // QE simulation
2962 //---
2963
2964 // set pointer to the QE table of the relevant pixel
2965
2966 qept = (float **)QE[ipixnum];
2967
2968 // check if wl is inside table; outside the table, QE is assumed to be zero
2969
2970 if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
2971 continue;
2972
2973 }
2974
2975 // find data point in the QE table (-> k)
2976
2977 k = 1; // start at 1 because the condition was already tested for 0
2978 while (k < pointsQE-1 && qept[0][k] < wl){
2979 k++;
2980 }
2981
2982 // calculate the qe value between 0. and 1.
2983
2984 qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
2985
2986 // if random > quantum efficiency, reject it
2987
2988 if ( (RandomNumber) > qe ) {
2989 continue;
2990 }
2991
2992 //+++
2993 // The photon has produced a photo electron
2994 //---
2995
2996 // cout << " accepted\n";
2997
2998 // increment the number of photoelectrons in the relevant pixel
2999
3000 nphe[ipixnum] += 1.0;
3001
3002 // store the new photoelectron
3003
3004 fadc->Fill(ipixnum,(t-*tmin_ns) , trigger->FillShow(ipixnum,t-*tmin_ns));
3005
3006 *itotnphe += 1;
3007 }
3008
3009 return(0);
3010
3011}
3012
3013
3014//------------------------------------------------------------
3015// @name produce_nsbrates
3016//
3017// @desc read the starfield file, call produce_phes on it in,
3018// @desc different wavebands, calculate the nsbrates
3019//
3020// @date Mon Feb 14 16:44:21 CET 2000
3021// @function @code
3022//------------------------------------------------------------
3023
3024int produce_nsbrates( char *iname, // the starfield input file name
3025 struct camera *cam, // camera layout
3026 float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
3027 // the NSB rates in phe/ns for each pixel
3028 ){
3029
3030 int i, j, k, ii; // counters
3031
3032 MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
3033 MFadc flashadc;
3034
3035 static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
3036 WAVEBANDBOUND2,
3037 WAVEBANDBOUND3,
3038 WAVEBANDBOUND4,
3039 WAVEBANDBOUND5,
3040 WAVEBANDBOUND6 };
3041
3042 FILE *infile; // the input file
3043 fpos_t fileposition; // marker on the input file
3044 static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
3045 static MCEventHeader evth; // the event header
3046 static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
3047 int itnphe; // total number of produced photoelectrons
3048 int itotnphe; // total number of produced photoelectrons after averaging
3049 int incph; // total number of cph read
3050 float tmin_ns; // minimum arrival time of all phes
3051 float tmax_ns; // maximum arrival time of all phes
3052 float integtime_ns; // integration time from the starfield generator
3053
3054 // open input file
3055
3056 log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
3057
3058 infile = fopen( iname, "r" );
3059
3060 // check if the satrfield input file exists
3061
3062 if ( infile == NULL ) {
3063
3064 log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
3065
3066 log( SIGNATURE, "There is not NSB from the Stars\n");
3067
3068 return (0);
3069 }
3070
3071 // get signature, and check it
3072
3073 if(check_reflector_file(infile)==FALSE){
3074 exit(1);
3075 }
3076
3077 // initialize flag
3078
3079 strcpy( cflag, " \0" );
3080
3081 // get flag
3082
3083 fread( cflag, SIZE_OF_FLAGS, 1, infile );
3084
3085 if ( ! feof(infile)){
3086
3087 // reading .rfl file
3088
3089 if(!isA( cflag, FLAG_START_OF_RUN )){
3090 error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
3091 }
3092 else { // found start of run
3093
3094 fread( cflag, SIZE_OF_FLAGS, 1, infile );
3095
3096 if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
3097
3098 // get MCEventHeader
3099
3100 fread( (char*)&evth, evth.mysize(), 1, infile );
3101
3102 integtime_ns = evth.get_energy();
3103
3104 // memorize where we are in the file
3105
3106 if (fgetpos( infile, &fileposition ) != 0){
3107 error( SIGNATURE, "Cannot position in file ...\n");
3108 }
3109
3110 // loop over the wavebands
3111
3112 for(i=0; i<iNUMWAVEBANDS; i++){
3113
3114 // initialize the rate array
3115
3116 for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
3117 rate_phepns[j][i] = 0.;
3118 }
3119
3120 itotnphe = 0;
3121
3122 // read the photons and produce the photoelectrons
3123 // - in order to average over the QE simulation, call the
3124 // production function iNUMNSBPRODCALLS times
3125
3126 for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
3127
3128 // position the file pointer to the beginning of the photons
3129
3130 fsetpos( infile, &fileposition);
3131
3132 // produce photoelectrons
3133
3134 k = produce_phes( infile,
3135 cam,
3136 wl_nm[i],
3137 wl_nm[i+1],
3138 &trigger, // this is a dummy here
3139 &flashadc, // this is a dummy here
3140 &itnphe,
3141 nphe, // we want this!
3142 &incph,
3143 &tmin_ns,
3144 &tmax_ns );
3145
3146 if( k != 0 ){ // non-zero returnvalue means error
3147 cout << "Exiting.\n";
3148 exit(1);
3149 }
3150
3151 for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
3152 rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
3153 }
3154
3155 itotnphe += itnphe;
3156
3157 } // end for(ii=0 ...
3158
3159 fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
3160 wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
3161
3162 } // end for(i=0 ...
3163
3164 }
3165 else{
3166 cout << "Starfield file contains no event.\nExiting.\n";
3167 exit(1);
3168 } // end if( isA ... event
3169 } // end if ( isA ... run
3170 }
3171 else{
3172 cout << "Starfield file contains no run.\nExiting.\n";
3173 exit(1);
3174 }
3175
3176 fclose( infile );
3177
3178 return(0);
3179
3180}
3181
3182
3183//------------------------------------------------------------
3184// @name produce_nsb_phes
3185//
3186// @desc produce the photoelectrons from the nsbrates
3187//
3188// @date Thu Feb 17 17:10:40 CET 2000
3189// @function @code
3190//------------------------------------------------------------
3191
3192int produce_nsb_phes( float *atmin_ns,
3193 float *atmax_ns,
3194 float theta_rad,
3195 struct camera *cam,
3196 float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
3197 float difnsbr_phepns[iMAXNUMPIX],
3198 float extinction[iNUMWAVEBANDS],
3199 float fnpx[iMAXNUMPIX],
3200 MTrigger *trigger,
3201 MFadc *fadc,
3202 int *inphe,
3203 float base_mv[iMAXNUMPIX]){
3204
3205 float simtime_ns; // NSB simulation time
3206 int i, j, k, ii;
3207 float zenfactor; // correction factor calculated from the extinction
3208 int inumnsbphe; // number of photoelectrons caused by NSB
3209
3210 float t;
3211 TRandom random; // Random numbers generator
3212
3213 random.SetSeed ((UInt_t) (RandomNumber*1000));
3214
3215 ii = *inphe; // avoid dereferencing
3216
3217 // check if the arrival times are set; if not generate them
3218
3219 if(*atmin_ns <SIMTIMEOFFSET_NS || *atmin_ns > *atmax_ns){
3220 *atmin_ns = 0.;
3221 *atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
3222
3223 }
3224 else{ // extend the event time window by the given offsets
3225
3226 *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
3227 *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
3228
3229 simtime_ns = *atmax_ns - *atmin_ns;
3230
3231 // make sure the simulated time is long enough for the FADC
3232 // simulation and not too long
3233
3234 if(simtime_ns< TOTAL_TRIGGER_TIME){
3235 *atmin_ns = *atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
3236 *atmax_ns = *atmin_ns + TOTAL_TRIGGER_TIME;
3237 simtime_ns = TOTAL_TRIGGER_TIME;
3238 }
3239
3240 if(simtime_ns> TOTAL_TRIGGER_TIME){
3241 *atmax_ns =*atmin_ns + TOTAL_TRIGGER_TIME;
3242 simtime_ns = TOTAL_TRIGGER_TIME;
3243 }
3244
3245 }
3246
3247 // initialize baselines
3248
3249 for(i=0; i<cam->inumpixels; i++){
3250 base_mv[i] = 0.;
3251 }
3252
3253 // calculate baselines and generate phes
3254
3255 for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
3256
3257 // calculate the effect of the atmospheric extinction
3258
3259 zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
3260
3261 for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
3262
3263 inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
3264 * simtime_ns );
3265
3266 base_mv[j] += inumnsbphe;
3267
3268 // randomize
3269
3270 if (inumnsbphe>0.0){
3271 inumnsbphe = random.Poisson (inumnsbphe );
3272 }
3273
3274 // create the photoelectrons
3275
3276 for(k=0; k < inumnsbphe; k++){
3277
3278 t=(RandomNumber * simtime_ns);
3279
3280 (*fadc).Fill(j,t ,(*trigger).FillNSB(j,t));
3281
3282 ii++; // increment total number of photoelectons
3283
3284 fnpx[j] += 1.; // increment number of photoelectrons in this pixel
3285
3286 }
3287
3288 } // end for(j=0 ...
3289 } // end for(i=0 ...
3290
3291 // finish baseline calculation
3292
3293 for(i=0; i<cam->inumpixels; i++){
3294 base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
3295 }
3296
3297 *inphe = ii; // update the pointer
3298
3299 return(0);
3300}
3301
3302
3303// @endcode
3304
3305
3306//=------------------------------------------------------------
3307//!@subsection Log of this file.
3308
3309//!@{
3310//
3311// $Log: not supported by cvs2svn $
3312// Revision 1.29 2001/11/14 17:38:23 blanch
3313// Sveral changes have been done:
3314// - bpoint_in_in_pixel has been dodified to speed up the program
3315// - Some minor changes have been done to remove warnings
3316// - buffer size and split version of the Branches have been removed
3317// - Some modifications were needed to adat the program to the new
3318// MRawEvtData::DeletePixels
3319//
3320// Revision 1.28 2001/10/26 16:31:45 blanch
3321// Removing several warnings.
3322//
3323// Revision 1.27 2001/09/05 10:04:33 blanch
3324// *** empty log message ***
3325//
3326// Revision 1.26 2001/07/19 09:29:53 blanch
3327// Different threshold for each pixel can be used.
3328//
3329// Revision 1.25 2001/05/08 08:07:54 blanch
3330// New numbering for branches from different trigger conditions has been
3331// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
3332// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
3333// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
3334// information already filled.
3335//
3336// Revision 1.24 2001/04/06 16:48:09 magicsol
3337// New camera version able to read the new format of the reflector output:
3338// reflector 0.4
3339//
3340// Revision 1.23 2001/03/28 16:13:41 blanch
3341// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
3342// been solved.
3343//
3344// Revision 1.22 2001/03/20 18:52:43 blanch
3345// *** empty log message ***
3346//
3347// Revision 1.21 2001/03/19 19:53:03 blanch
3348// Some print outs have been removed.
3349//
3350// Revision 1.20 2001/03/19 19:30:06 magicsol
3351// Minor changes have been done to improve the FADC pedestals treatment.
3352//
3353// Revision 1.19 2001/03/05 11:14:41 magicsol
3354// I changed the position of readinf a parameter. It is a minnor change.
3355//
3356// Revision 1.18 2001/03/05 10:36:52 blanch
3357// A branch with information about the FADC simulation (MMcFadcHeader) is writen
3358// in the McHeaders tree of the aoutput root file.
3359// The NSB contribution is only applied if the the number of phe form the shower
3360// are above a given number.
3361//
3362// Revision 1.17 2001/02/23 11:05:57 magicsol
3363// Small changes due to slightly different output format and the introduction of
3364// pedesals for teh FADC.
3365//
3366// Revision 1.16 2001/01/18 18:44:40 magicsol
3367// *** empty log message ***
3368//
3369// Revision 1.15 2001/01/17 09:32:27 harald
3370// The changes are neccessary to have the same name for trees in MC and in
3371// data. So now there should be no differences in MC data and real data from
3372// FADC system.
3373//
3374// Revision 1.14 2001/01/15 12:33:34 magicsol
3375// Some modifications have been done to use the new (Dec'2000) Raw data format.
3376// There are also some minnors modifications to adapt the improvements in the
3377// MTrigger class (overlaping time and trigger cells).
3378//
3379// Revision 1.13 2000/10/25 08:14:23 magicsol
3380// The routine that produce poisson random numbers to decide how many phe
3381// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
3382// is used.
3383//
3384// Revision 1.12 2000/09/22 17:40:18 harald
3385// Added a lot of changes done by oscar.
3386//
3387// Revision 1.11 2000/09/21 11:47:33 harald
3388// Oscar found some smaller errors in the calculation of the pixel shape and
3389// corrected it.
3390//
3391// $Log: not supported by cvs2svn $
3392// Revision 1.29 2001/11/14 17:38:23 blanch
3393// Sveral changes have been done:
3394// - bpoint_in_in_pixel has been dodified to speed up the program
3395// - Some minor changes have been done to remove warnings
3396// - buffer size and split version of the Branches have been removed
3397// - Some modifications were needed to adat the program to the new
3398// MRawEvtData::DeletePixels
3399//
3400// Revision 1.28 2001/10/26 16:31:45 blanch
3401// Removing several warnings.
3402//
3403// Revision 1.27 2001/09/05 10:04:33 blanch
3404// *** empty log message ***
3405//
3406// Revision 1.26 2001/07/19 09:29:53 blanch
3407// Different threshold for each pixel can be used.
3408//
3409// Revision 1.25 2001/05/08 08:07:54 blanch
3410// New numbering for branches from different trigger conditions has been
3411// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
3412// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
3413// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
3414// information already filled.
3415//
3416// Revision 1.24 2001/04/06 16:48:09 magicsol
3417// New camera version able to read the new format of the reflector output:
3418// reflector 0.4
3419//
3420// Revision 1.23 2001/03/28 16:13:41 blanch
3421// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
3422// been solved.
3423//
3424// Revision 1.22 2001/03/20 18:52:43 blanch
3425// *** empty log message ***
3426//
3427// Revision 1.21 2001/03/19 19:53:03 blanch
3428// Some print outs have been removed.
3429//
3430// Revision 1.20 2001/03/19 19:30:06 magicsol
3431// Minor changes have been done to improve the FADC pedestals treatment.
3432//
3433// Revision 1.19 2001/03/05 11:14:41 magicsol
3434// I changed the position of readinf a parameter. It is a minnor change.
3435//
3436// Revision 1.18 2001/03/05 10:36:52 blanch
3437// A branch with information about the FADC simulation (MMcFadcHeader) is writen
3438// in the McHeaders tree of the aoutput root file.
3439// The NSB contribution is only applied if the the number of phe form the shower
3440// are above a given number.
3441//
3442// Revision 1.17 2001/02/23 11:05:57 magicsol
3443// Small changes due to slightly different output format and the introduction of
3444// pedesals for teh FADC.
3445//
3446// Revision 1.16 2001/01/18 18:44:40 magicsol
3447// *** empty log message ***
3448//
3449// Revision 1.15 2001/01/17 09:32:27 harald
3450// The changes are neccessary to have the same name for trees in MC and in
3451// data. So now there should be no differences in MC data and real data from
3452// FADC system.
3453//
3454// Revision 1.14 2001/01/15 12:33:34 magicsol
3455// Some modifications have been done to use the new (Dec'2000) Raw data format.
3456// There are also some minnors modifications to adapt the improvements in the
3457// MTrigger class (overlaping time and trigger cells).
3458//
3459// Revision 1.13 2000/10/25 08:14:23 magicsol
3460// The routine that produce poisson random numbers to decide how many phe
3461// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
3462// is used.
3463//
3464// Revision 1.10 2000/07/04 14:10:20 MagicSol
3465// Some changes have been done in the root output file. The RawEvt tree is only
3466// stored in the single trigger mode.
3467// The trigger input parameters are also given by the general input card.
3468// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
3469// each one can be studied seperately.
3470//
3471// Revision 1.9 2000/06/13 13:25:24 blanch
3472// The multiple arrays have been replaced, since they do not work
3473// in alpha machines. Now we are using pointers and the command new
3474// to allocate memory.
3475//
3476// Revision 1.8 2000/05/11 13:57:27 blanch
3477// The option to loop over several trigger configurations has been included.
3478// Some non-sense with triggertime range has been solved.
3479// Montecarlo information and ADC counts are saved in a root file.
3480// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
3481//
3482// Revision 1.7 2000/03/24 18:10:46 blanch
3483// A first FADC simulation and a trigger simulation are already implemented.
3484// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
3485// A loop over trigger threshold and some corretcions in the time range where it looks for a trigger will be implemented as soon as possible.
3486//
3487// Revision 1.6 2000/03/20 18:35:11 blanch
3488// The trigger is already implemented but it does not save the trigger information in any file as it is implemented in timecam. In the next days there will be a version which also creates the files with the trigger information. It is going to be a mixing of the current camera and timecam programs.
3489//
3490// Revision 1.5 2000/02/18 17:40:35 petry
3491// This version includes drastic changes compared to camera.cxx 1.4.
3492// It is not yet finished and not immediately useful because the
3493// trigger simulation is not yet re-implemented. I had to take it
3494// out together with some other stuff in order to tidy the whole
3495// program up. This is not meant as an insult to anyone. I needed
3496// to do this in order to be able to work on it.
3497//
3498// This version has been put in the repository in order to be
3499// able to share the further development with others.
3500//
3501// If you need something working, wait or take an earlier one.
3502// See file README.
3503//
3504// Revision 1.4 2000/01/25 08:36:23 petry
3505// The pixelization in previous versions was buggy.
3506// This is the first version with a correct pixelization.
3507//
3508// Revision 1.3 2000/01/20 18:22:17 petry
3509// Found little bug which makes camera crash if it finds a photon
3510// of invalid wavelength. This bug is now fixed and the range
3511// of valid wavelengths extended to 290 - 800 nm.
3512// This is in preparation for the NSB simulation to come.
3513// Dirk
3514//
3515// Revision 1.2 1999/11/19 08:40:42 harald
3516// Now it is possible to compile the camera programm under osf1.
3517//
3518// Revision 1.1.1.1 1999/11/05 11:59:31 harald
3519// This the starting point for CVS controlled further developments of the
3520// camera program. The program was originally written by Jose Carlos.
3521// But here you can find a "rootified" version to the program. This means
3522// that there is no hbook stuff in it now. Also the output of the
3523// program changed to the MagicRawDataFormat.
3524//
3525// The "rootification" was done by Dirk Petry and Harald Kornmayer.
3526//
3527// Revision 1.3 1999/10/22 15:01:28 petry
3528// version sent to H.K. and N.M. on Fri Oct 22 1999
3529//
3530// Revision 1.2 1999/10/22 09:44:23 petry
3531// first synthesized version which compiles and runs without crashing;
3532//
3533// Revision 1.1.1.1 1999/10/21 16:35:10 petry
3534// first synthesised version
3535//
3536// Revision 1.13 1999/03/15 14:59:05 gonzalez
3537// camera-1_1
3538//
3539// Revision 1.12 1999/03/02 09:56:10 gonzalez
3540// *** empty log message ***
3541//
3542//
3543//!@}
3544
3545//=EOF
Note: See TracBrowser for help on using the repository browser.