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

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