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

Last change on this file since 1598 was 1597, checked in by blanch, 23 years ago
*** empty log message ***
File size: 111.6 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.46 $
24// $Author: blanch $
25// $Date: 2002-11-08 17:51:00 $
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 Float_t ftime, ltime;
1591 if (reflector_file_version<6){
1592 mcevth.get_times(&ftime, &ltime);
1593 McEvt->Fill( -1.0,
1594 (UShort_t) mcevth.get_primary() ,
1595 mcevth.get_energy(),
1596 -1.0,
1597 -1.0,
1598 -1.0,
1599 mcevth.get_theta(),
1600 mcevth.get_phi(),
1601 mcevth.get_core(),
1602 mcevth.get_coreX(),
1603 mcevth.get_coreY(),
1604 impactD,
1605 mcevth_2.get_theta_CT(),
1606 mcevth_2.get_phi_CT(),
1607 ftime,
1608 ltime,
1609 0,
1610 0,
1611 0,
1612 0,
1613 0,
1614 0,
1615 0,
1616 (UInt_t)mcevth.get_CORSIKA(),
1617 (UInt_t)mcevth.get_AtmAbs(),
1618 (UInt_t)mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot(),
1619 (UInt_t) ncph,
1620 (UInt_t) inumphe,
1621 (UInt_t) inumphensb+inumphe);
1622 }
1623 else{
1624 Float_t Nmax, t0, tmax, a, b, c, chi2;
1625 mcevth.get_times(&ftime, &ltime);
1626 chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
1627 McEvt->Fill( mcevth_2.get_evt_number(),
1628 (UShort_t) mcevth_2.get_primary() ,
1629 mcevth_2.get_energy(),
1630 mcevth_2.get_thick0(),
1631 mcevth_2.get_first_target(),
1632 mcevth_2.get_z_first_int(),
1633 mcevth_2.get_theta(),
1634 mcevth_2.get_phi(),
1635 mcevth_2.get_core(),
1636 mcevth_2.get_coreX(),
1637 mcevth_2.get_coreY(),
1638 impactD,
1639 mcevth_2.get_theta_CT(),
1640 mcevth_2.get_phi_CT(),
1641 ftime,
1642 ltime,
1643 Nmax,
1644 t0,
1645 tmax,
1646 a,
1647 b,
1648 c,
1649 chi2,
1650 (UInt_t)mcevth_2.get_CORSIKA(),
1651 (UInt_t)mcevth_2.get_AtmAbs(),
1652 (UInt_t)mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot(),
1653 (UInt_t) ncph,
1654 (UInt_t) inumphe,
1655 (UInt_t) inumphensb+inumphe);
1656 }
1657 }
1658 // Fill the Tree with the current leaves of each branch
1659 i=EvtTree.Fill() ;
1660
1661 // Clear the branches
1662 if(Write_McTrig){
1663 for(i=0;i<icontrigger;i++){
1664 McTrig[i]->Clear() ;
1665 }
1666 }
1667 if( Write_RawEvt ){
1668 for(i=0;i<icontrigger;i++){
1669 EvtHeader[i]->Clear() ;
1670 EvtData[i]->DeletePixels();
1671 }
1672 }
1673 if (Write_McEvt)
1674 McEvt->Clear() ;
1675 }
1676 }
1677 // We study a single trigger condition
1678 else {
1679
1680 // Setting trigger conditions
1681 Trigger.SetMultiplicity(Trigger_multiplicity);
1682 Trigger.SetTopology(Trigger_topology);
1683 for (int i=0;i<CAMERA_PIXELS;i++)
1684 fpixelthres[i]=qThreshold[i];
1685 Trigger.SetThreshold(fpixelthres);
1686
1687 Trigger.Diskriminate() ;
1688
1689 //
1690 // look if in all the signals in the trigger signal branch
1691 // is a possible Trigger. Therefore we habe to diskriminate all
1692 // the simulated analog signals (Method Diskriminate in class
1693 // MTrigger). We look simultanously for the moments at which
1694 // there are more than TRIGGER_MULTI pixels above the
1695 // CHANNEL_THRESHOLD.
1696 //
1697
1698 Lev0 = (Short_t) Trigger.ZeroLevel() ;
1699
1700 Lev1 = 0 ;
1701
1702 //
1703 // Start the First Level Trigger simulation
1704 //
1705
1706 if ( Lev0 > 0 || Write_All_Images) {
1707 Lev1= Trigger.FirstLevel();
1708 }
1709 if (Lev1>0){
1710 ++ntrigger;
1711 }
1712 Int_t NumImages = Lev1;
1713 Lev0=1;
1714 if (Lev1==0 && Write_All_Images){
1715 NumImages=1;
1716 Lev0=0;
1717 }
1718
1719 for(Int_t ii=0;ii<NumImages;ii++){
1720 // Loop over different level one triggers
1721
1722 //
1723 // fill inside class fadc the member output
1724 //
1725 fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
1726
1727 if (Write_McTrig){
1728 McTrig[0]->SetFirstLevel ((ii+1)*Lev0);
1729 McTrig[0]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
1730 Trigger.GetMapDiskriminator(trigger_map);
1731 McTrig[0]->SetMapPixels(trigger_map,ii);
1732 }
1733
1734 // Fill Evt information
1735
1736 if (Write_RawEvt){
1737
1738 //
1739 // Fill the header of this event
1740 //
1741
1742 EvtHeader[0]->FillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ;
1743
1744 // fill pixel information
1745
1746 if (Lev1){
1747 for(UInt_t i=0;i<camgeom.GetNumPixels();i++){
1748 for (j=0;j<FADC_SLICES;j++){
1749 fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
1750 }
1751 EvtData[0]->AddPixel(i,fadcValues,0);
1752 }
1753 }
1754 }
1755 //
1756 // fill the MMcEvt with all information
1757 //
1758
1759 if (Write_McEvt){
1760 }
1761 // We don not count photons out of the camera.
1762
1763
1764 //
1765 // write it out to the file outfile
1766 //
1767
1768 EvtTree.Fill() ;
1769
1770 //
1771 // if a first level trigger occurred, then
1772 // 1. do some other stuff (not implemented)
1773 // 2. start the gui tool
1774
1775 if(FADC_Scan){
1776 if ( Lev0 > 0 ) {
1777 fadc.ShowSignal( McEvt, (Float_t) 60. ) ;
1778 }
1779 }
1780
1781 if(Trigger_Scan){
1782 if ( Lev0 > 0 ) {
1783 Trigger.ShowSignal(McEvt) ;
1784 }
1785 }
1786
1787 // clear all
1788 if (Write_RawEvt) EvtHeader[0]->Clear() ;
1789 if (Write_RawEvt) EvtData[0]->DeletePixels();
1790 if (Write_McEvt) McEvt->Clear() ;
1791 }
1792 if (Write_McTrig) McTrig[0]->Clear() ;
1793 }
1794
1795#ifdef __DEBUG__
1796 printf("\n");
1797
1798 for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1799
1800 for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1801
1802 if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1803
1804 if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1805
1806 printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1807 (int)pixels[ici][icj][PIXNUM],
1808 pixels[ici][icj][PIXX],
1809 pixels[ici][icj][PIXY],
1810 fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1811
1812 }
1813 }
1814 }
1815 }
1816
1817 for (int i=0; i<camgeom.GetNumPixels(); ++i) {
1818 printf("%d (%d): ", i, npixneig[i]);
1819 for (j=0; j<npixneig[i]; ++i)
1820 printf(" %d", pixneig[i][j]);
1821 printf("\n");
1822 }
1823
1824#endif // __DEBUG__
1825
1826
1827 // look for the next event
1828
1829 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1830
1831 } // end while there is a next event
1832
1833 if( !isA( flag, FLAG_END_OF_RUN )){
1834 error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
1835 }
1836 else { // found end of run
1837 ntshow += nshow;
1838 log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
1839
1840 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1841
1842 if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
1843 log(SIGNATURE, "End of file . . .\n");
1844 still_in_loop = FALSE;
1845
1846 if ((! Data_From_STDIN) && ( !feof(inputfile) )){
1847
1848 // we have concatenated input files.
1849 // get signature of the next part and check it.
1850
1851 if((reflector_file_version=check_reflector_file( inputfile ))==FALSE){
1852 exit(1);
1853 }
1854
1855 }
1856
1857 fread( flag, SIZE_OF_FLAGS, 1, inputfile );
1858
1859 } // end if found end of file
1860 } // end if found end of run
1861
1862 } // end if else found start of run
1863 } // end big while loop
1864
1865 //<@ Finally we should fill th McRunHeader
1866
1867 Float_t heights[10];
1868 time_t ltime;
1869 Float_t ftime;
1870 Float_t rnum;
1871
1872 get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
1873 get_teles_axis(&telestheta,&telesphi);
1874 get_source_off(&sofftheta,&soffphi);
1875 mcevth.get_theta_range(&shthetamin, &shthetamax);
1876 mcevth.get_phi_range(&shphimin,&shphimax);
1877 corsika=UInt_t(mcevth.get_VersionPGM()*1000);
1878 if(!Trigger_Loop) icontrigger=0;
1879 time (&ltime);
1880 ftime = ((Float_t)ltime)/1000;
1881 for (int i=0; i< 10;i++)
1882 heights[i]=mcevth.get_HeightLev (i);
1883 rnum=mcevth.get_RunNumber();
1884
1885 McRunHeader->Fill(rnum,
1886 (UInt_t) 0,
1887 mcevth.get_DateRun(),
1888 ftime,
1889 icontrigger,
1890 !Write_All_Images,
1891 Write_McEvt,
1892 Write_McTrig,
1893 Write_McFADC,
1894 Write_RawEvt,
1895 addElecNoise,
1896 CAMERA_PIXELS,
1897 (UInt_t)ntshow,
1898 (UInt_t)ntrigger,
1899 0,
1900 sfRaH,
1901 sfRaM,
1902 sfRaS,
1903 sfDeD,
1904 sfDeM,
1905 sfDeS,
1906 meanNSB,
1907 telestheta,
1908 telesphi,
1909 sofftheta,
1910 soffphi,
1911 shthetamax,
1912 shthetamin,
1913 shphimax,
1914 shphimin,
1915 mcevth.get_CWaveLower(),
1916 mcevth.get_CWaveUpper(),
1917 mcevth.get_slope(),
1918 mcevth.get_NumObsLev(),
1919 heights,
1920 corsika,
1921 (UInt_t)(REFL_VERSION_A*100),
1922 (UInt_t)(VERSION*100),
1923 0);
1924
1925 // Fill the Header Tree with the current leaves of each branch
1926 HeaderTree.Fill() ;
1927
1928 //++
1929 // put the Event to the root file
1930 //--
1931
1932 outfile_temp.Write() ;
1933 outfile_temp.Close() ;
1934
1935 // close input file
1936
1937 log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
1938 ntshow, ntcph );
1939 datafile<<ntshow<<" event(s), with a total of "<<ntcph<<" C.photons"<<endl;
1940 if (Trigger_Loop){
1941 log( SIGNATURE, "Trigger Mode. \n");
1942 log( SIGNATURE, "Fraction of triggers: \n");
1943 datafile<<"Fraction of triggers: "<<endl;
1944 for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1945 for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1946 for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1947 log( SIGNATURE, "Thres %5.1f, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
1948 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);
1949 datafile<<"Thres "<<fthrescount<<", Multi "<<imulticount<<", Topo"<<isorttopo[itopocount]<<": ";
1950 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;
1951 }
1952 }
1953 }
1954 }
1955 else{
1956 log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
1957 ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
1958 datafile<<"Fraction of triggers: "<<((float)ntrigger) / ((float)ntshow) * 100.0<<" ("<<ntrigger<<" out of "<<ntshow<<" )"<<endl;
1959 }
1960
1961 // close files
1962
1963 log( SIGNATURE, "Closing files\n" );
1964
1965 if( ! Data_From_STDIN ){
1966 fclose( inputfile );
1967 }
1968 datafile.close();
1969
1970 // program finished
1971
1972 log( SIGNATURE, "Done.\n");
1973
1974 return( 0 );
1975}
1976//!@}
1977
1978// @T \newpage
1979
1980//!@subsection Functions definition.
1981
1982//!-----------------------------------------------------------
1983// @name present
1984//
1985// @desc Make some presentation
1986//
1987// @date Sat Jun 27 05:58:56 MET DST 1998
1988//------------------------------------------------------------
1989// @function
1990
1991//!@{
1992void
1993present(void)
1994{
1995 cout << "##################################################\n"
1996 << SIGNATURE << '\n' << '\n'
1997 << "Processor of the reflector output\n"
1998 << "J C Gonzalez, Jun 1998\n"
1999 << "##################################################\n\n"
2000 << flush ;
2001}
2002//!@}
2003
2004
2005//!-----------------------------------------------------------
2006// @name usage
2007//
2008// @desc show help
2009//
2010// @date Tue Dec 15 16:23:30 MET 1998
2011//------------------------------------------------------------
2012// @function
2013
2014//!@{
2015void
2016usage(void)
2017{
2018 present();
2019 cout << "\nusage ::\n\n"
2020 << "\t camera "
2021 << " [ -@ paramfile ] "
2022 << " [ -h ] "
2023 << "\n\n or \n\n"
2024 << "\t camera < paramfile"
2025 << "\n\n";
2026 exit(0);
2027}
2028//!@}
2029
2030
2031//!-----------------------------------------------------------
2032// @name log
2033//
2034// @desc function to send log information
2035//
2036// @var funct Name of the caller function
2037// @var fmt Format to be used (message)
2038// @var ... Other information to be shown
2039//
2040// @date Sat Jun 27 05:58:56 MET DST 1998
2041//------------------------------------------------------------
2042// @function
2043
2044//!@{
2045void
2046log(const char *funct, char *fmt, ...)
2047{
2048 va_list args;
2049
2050 // Display the name of the function that called error
2051 printf("[%s]: ", funct);
2052
2053 // Display the remainder of the message
2054 va_start(args, fmt);
2055 vprintf(fmt, args);
2056 va_end(args);
2057}
2058//!@}
2059
2060
2061//!-----------------------------------------------------------
2062// @name error
2063//
2064// @desc function to send an error message, and abort the program
2065//
2066// @var funct Name of the caller function
2067// @var fmt Format to be used (message)
2068// @var ... Other information to be shown
2069//
2070// @date Sat Jun 27 05:58:56 MET DST 1998
2071//------------------------------------------------------------
2072// @function
2073
2074//!@{
2075void
2076error(const char *funct, char *fmt, ...)
2077{
2078 va_list args;
2079
2080 // Display the name of the function that called error
2081 fprintf(stdout, "ERROR in %s: ", funct);
2082
2083 // Display the remainder of the message
2084 va_start(args, fmt);
2085 vfprintf(stdout, fmt, args);
2086 va_end(args);
2087
2088 perror(funct);
2089
2090 exit(1);
2091}
2092//!@}
2093
2094
2095//!-----------------------------------------------------------
2096// @name isA
2097//
2098// @desc returns TRUE(FALSE), if the flag is(is not) the given
2099//
2100// @var s1 String to be searched
2101// @var flag Flag to compare with string s1
2102// @return TRUE: both strings match; FALSE: oth.
2103//
2104// @date Wed Jul 8 15:25:39 MET DST 1998
2105//------------------------------------------------------------
2106// @function
2107
2108//!@{
2109int
2110isA( char * s1, const char * flag ) {
2111 return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
2112}
2113//!@}
2114
2115
2116//!-----------------------------------------------------------
2117// @name read_ct_file
2118//
2119// @desc read CT definition file
2120//
2121// @date Sat Jun 27 05:58:56 MET DST 1998
2122//------------------------------------------------------------
2123// @function
2124
2125//!@{
2126void
2127read_ct_file(void)
2128{
2129 char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
2130 char token[ITEM_MAX_LENGTH]; //@< a single token
2131 int i, j; //@< dummy counters
2132
2133 log( "read_ct_file", "start.\n" );
2134
2135 ifstream ctin ( ct_filename );
2136
2137 if ( ctin.bad() )
2138 error( "read_ct_file",
2139 "Cannot open CT def. file: %s\n", ct_filename );
2140
2141 // loop till the "end" directive is reached
2142
2143 while (!ctin.eof()) {
2144
2145 // get line from stdin
2146
2147 ctin.getline(line, LINE_MAX_LENGTH);
2148
2149 // look for each item at the beginning of the line
2150
2151 for (i=0; i<=define_mirrors; i++)
2152 if (strstr(line, CT_ITEM_NAMES[i]) == line)
2153 break;
2154
2155 // if it is not a valid line, just ignore it
2156
2157 if (i == define_mirrors+1)
2158 continue;
2159
2160 // case block for each directive
2161
2162 switch ( i ) {
2163
2164 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2165
2166 // get focal distance
2167
2168 sscanf(line, "%s %d", token, &ct_Type);
2169
2170 log( "read_ct_file", "<Type of Telescope>: %s\n",
2171 ((ct_Type==0) ? "CT1" : "MAGIC") );
2172
2173 break;
2174
2175 case focal_distance: // <focal distance> [cm]
2176
2177 // get focal distance
2178
2179 sscanf(line, "%s %f", token, &ct_Focal_mean);
2180
2181 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2182
2183 break;
2184
2185 case focal_std: // s(focal distance) [cm]
2186
2187 // get focal distance
2188
2189 sscanf(line, "%s %f", token, &ct_Focal_std);
2190
2191 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2192
2193 break;
2194
2195 case point_spread: // <point spread> [cm]
2196
2197 // get point spread
2198
2199 sscanf(line, "%s %f", token, &ct_PSpread_mean);
2200
2201 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2202
2203 break;
2204
2205 case point_std: // s(point spread) [cm]
2206
2207 // get point spread
2208
2209 sscanf(line, "%s %f", token, &ct_PSpread_std);
2210
2211 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2212
2213 break;
2214
2215 case adjustment_dev: // s(adjustment_dev) [cm]
2216
2217 // get point spread
2218
2219 sscanf(line, "%s %f", token, &ct_Adjustment_std);
2220
2221 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2222
2223 break;
2224
2225 case black_spot: // radius of the black spot in the center [cm]
2226
2227 // get black spot radius
2228
2229 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2230
2231 log( "read_ct_file", "Radius of the black spots: %f cm\n",
2232 ct_BlackSpot_rad);
2233
2234 break;
2235
2236 case r_mirror: // radius of the mirrors [cm]
2237
2238 // get radius of mirror
2239
2240 sscanf(line, "%s %f", token, &ct_RMirror);
2241
2242 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2243
2244 break;
2245
2246 case n_mirrors: // number of mirrors
2247
2248 // get the name of the output_file from the line
2249
2250 sscanf(line, "%s %d", token, &ct_NMirrors);
2251
2252 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2253
2254 break;
2255
2256 case camera_width: // camera width [cm]
2257
2258 // get the name of the ct_file from the line
2259
2260 sscanf(line, "%s %f", token, &ct_CameraWidth);
2261
2262 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2263
2264 break;
2265
2266 case n_pixels: // number of pixels
2267
2268 // get the name of the output_file from the line
2269
2270 sscanf(line, "%s %d", token, &ct_NPixels);
2271
2272 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2273
2274 break;
2275
2276 case n_centralpixels: // number of central pixels
2277
2278 // get the name of the output_file from the line
2279
2280 sscanf(line, "%s %d", token, &ct_NCentralPixels);
2281
2282 log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
2283
2284 break;
2285
2286 case n_gappixels: // number of gap pixels
2287
2288 // get the name of the output_file from the line
2289
2290 sscanf(line, "%s %d", token, &ct_NGapPixels);
2291
2292 log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
2293
2294 break;
2295
2296 case pixel_width: // pixel width [cm]
2297
2298 // get the name of the ct_file from the line
2299
2300 sscanf(line, "%s %f", token, &ct_PixelWidth);
2301
2302 ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
2303
2304 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2305
2306 break;
2307
2308 case define_mirrors: // read table with the parameters of the mirrors
2309
2310 log( "read_ct_file", "Table of mirrors data:\n" );
2311
2312 // check whether the number of mirrors was already set
2313
2314 if ( ct_NMirrors == 0 )
2315 error( "read_ct_file", "NMirrors was not set.\n" );
2316
2317 // allocate memory for paths list
2318
2319 log( "read_ct_file", "Allocating memory for ct_data\n" );
2320
2321 ct_data = new float*[ct_NMirrors];
2322
2323 for (i=0; i<ct_NMirrors; i++)
2324 ct_data[i] = new float[CT_NDATA];
2325
2326 // read data
2327
2328 log( "read_ct_file", "Reading mirrors data...\n" );
2329
2330 for (i=0; i<ct_NMirrors; i++)
2331 for (j=0; j<CT_NDATA; j++)
2332 ctin >> ct_data[i][j];
2333
2334 break;
2335
2336 } // switch ( i )
2337
2338 } // end while
2339
2340 // end
2341
2342 log( "read_ct_file", "done.\n" );
2343
2344 return;
2345}
2346//!@}
2347
2348
2349//!-----------------------------------------------------------
2350// @name read_QE
2351//
2352// @desc read QE data
2353//
2354// @date thu 5 17:59:57 CEST 2002
2355//------------------------------------------------------------
2356// @function
2357
2358//!@{
2359void
2360read_QE(void){
2361 ifstream qefile;
2362 char line[LINE_MAX_LENGTH];
2363 int i, j, icount;
2364 float qe;
2365
2366 //------------------------------------------------------------
2367 // second, pixels' QE
2368
2369 // try to open the file
2370
2371 log("read_QE", "Opening the file \"%s\" . . .\n", QE_FILE);
2372
2373 qefile.open( QE_FILE );
2374
2375 // if it is wrong or does not exist, exit
2376
2377 if ( qefile.bad() )
2378 error( "read_QE", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2379
2380 // read file
2381
2382 log("read_QE", "Reading data . . .\n");
2383
2384 i=-1;
2385 icount = 0;
2386
2387 while ( ! qefile.eof() ) {
2388
2389 // get line from the file
2390
2391 qefile.getline(line, LINE_MAX_LENGTH);
2392
2393 // skip if comment
2394
2395 if ( *line == '#' )
2396 continue;
2397
2398 // if it is the first valid value, it is the number of QE data points
2399
2400 if ( i < 0 ) {
2401
2402 // get the number of datapoints
2403
2404 sscanf(line, "%d", &pointsQE);
2405
2406 // allocate memory for the table of QEs
2407
2408 QE = new float ** [ct_NPixels];
2409
2410 for ( i=0; i<ct_NPixels; ++i ) {
2411 QE[i] = new float * [2];
2412 QE[i][0] = new float[pointsQE];
2413 QE[i][1] = new float[pointsQE];
2414 }
2415
2416 QElambda = new float [pointsQE];
2417
2418 for ( i=0; i<pointsQE; ++i ) {
2419 qefile.getline(line, LINE_MAX_LENGTH);
2420 sscanf(line, "%f", &QElambda[i]);
2421 }
2422
2423 i=0;
2424
2425 continue;
2426 }
2427
2428 // get the values (num-pixel, num-datapoint, QE-value)
2429
2430 if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2431 break;
2432
2433 if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2434 ((j-1) < pointsQE) && ((j-1) > -1) ) {
2435 QE[i-1][0][j-1] = QElambda[j-1];
2436 QE[i-1][1][j-1] = qe;
2437 }
2438
2439 if ( i > ct_NPixels) break;
2440
2441 icount++;
2442
2443 }
2444
2445 if(icount/pointsQE < ct_NPixels){
2446 error( "read_QE", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2447 icount/pointsQE, ct_NPixels );
2448 }
2449
2450 // close file
2451
2452 qefile.close();
2453
2454 // test QE
2455
2456 for(icount=0; icount< ct_NPixels; icount++){
2457 for(i=0; i<pointsQE; i++){
2458 if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
2459 QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
2460 error( "read_QE", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2461 icount, i, QE[icount][0][i], QE[icount][1][i] );
2462 }
2463 }
2464 }
2465
2466 // end
2467
2468 log("read_QE", "Done.\n");
2469}
2470//!@}
2471
2472//!-----------------------------------------------------------
2473// @name read_WC
2474//
2475// @desc read WC data
2476//
2477// @date thu 5 17:59:57 CEST 2002
2478//------------------------------------------------------------
2479// @function
2480
2481//!@{
2482void
2483read_WC(void){
2484 ifstream wcfile;
2485 char line[LINE_MAX_LENGTH];
2486 int i;
2487
2488 //------------------------------------------------------------
2489 // Read Light Guides data
2490
2491 // try to open the file
2492
2493 log("read_QE", "Opening the file \"%s\" . . .\n", WC_FILE);
2494
2495 wcfile.open( WC_FILE );
2496
2497 // if it is wrong or does not exist, exit
2498
2499 if ( wcfile.bad() )
2500 error( "read_WC", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2501
2502 // read file
2503
2504 log("read_WC", "Reading data . . .\n");
2505
2506 // get line from the file
2507
2508 wcfile.getline(line, LINE_MAX_LENGTH);
2509
2510 // get the number of datapoints
2511
2512 sscanf(line, "%d", &pointsWC);
2513
2514 // allocate memory for the table of QEs
2515
2516 WC = new float * [2];
2517 WC[0] = new float[pointsWC];
2518 WC[1] = new float[pointsWC];
2519
2520 for ( i=0; i<pointsWC; ++i ) {
2521 wcfile.getline(line, LINE_MAX_LENGTH);
2522 sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
2523 }
2524
2525 // close file
2526
2527 wcfile.close();
2528
2529 // read
2530
2531 log("read_WC", "Done.\n");
2532}
2533//!@}
2534
2535//!-----------------------------------------------------------
2536// @name read_pixels
2537//
2538// @desc read pixels data
2539//
2540// @date Fri Mar 12 16:33:34 MET 1999
2541//------------------------------------------------------------
2542// @function
2543
2544//!@{
2545void
2546read_pixels(struct camera *pcam)
2547{
2548 ifstream qefile;
2549 char line[LINE_MAX_LENGTH];
2550 int n, i, j, icount;
2551 float qe;
2552
2553 //------------------------------------------------------------
2554 // first, pixels' coordinates
2555
2556 pcam->inumpixels = ct_NPixels;
2557 pcam->inumcentralpixels = ct_NCentralPixels;
2558 pcam->inumgappixels = ct_NGapPixels;
2559 pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
2560 pcam->dpixdiameter_cm = ct_PixelWidth;
2561
2562 // initialize pixel numbers
2563
2564 pixary = new float* [2*ct_NCentralPixels];
2565 for ( i=0; i<2*ct_NCentralPixels; ++i )
2566 pixary[i] = new float[2];
2567
2568 pixneig = new int* [ct_NCentralPixels];
2569 for ( i=0; i<ct_NCentralPixels; ++i ) {
2570 pixneig[i] = new int[6];
2571 for ( j=0; j<6; ++j )
2572 pixneig[i][j] = -1;
2573 }
2574
2575 npixneig = new int[ct_NCentralPixels];
2576 for ( i=0; i<ct_NCentralPixels; ++i )
2577 npixneig[i] = 0;
2578
2579 // generate all coordinates
2580
2581 igen_pixel_coordinates(pcam);
2582
2583
2584 // calculate tables of neighbours
2585
2586#ifdef __DEBUG__
2587 for ( n=0 ; n<ct_NPixels ; ++n ) {
2588 cout << "Para el pixel " << n << ": ";
2589 for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2590 if ( pixels_are_neig(n,i) == TRUE ) {
2591 pixneig[n][npixneig[n]] = i;
2592 pixneig[i][npixneig[i]] = n;
2593 cout << i << ' ';
2594 ++npixneig[n];
2595 ++npixneig[i];
2596 }
2597 }
2598 cout << endl << flush;
2599 }
2600#else // ! __DEBUG__
2601 for ( n=0 ; n<ct_NCentralPixels ; ++n )
2602 for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2603 if ( pixels_are_neig(n,i) == TRUE ) {
2604 pixneig[n][npixneig[n]] = i;
2605 pixneig[i][npixneig[i]] = n;
2606 ++npixneig[n];
2607 ++npixneig[i];
2608 }
2609#endif // ! __DEBUG__
2610
2611#ifdef __DEBUG__
2612 for ( n=0 ; n<ct_NPixels ; ++n ) {
2613 cout << n << ':';
2614 for ( j=0; j<npixneig[n]; ++j)
2615 cout << ' ' << pixneig[n][j];
2616 cout << endl << flush;
2617 }
2618#endif // __DEBUG__
2619
2620 //------------------------------------------------------------
2621 // second, pixels' QE
2622
2623 // try to open the file
2624
2625 log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
2626
2627 qefile.open( QE_FILE );
2628
2629 // if it is wrong or does not exist, exit
2630
2631 if ( qefile.bad() )
2632 error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2633
2634 // read file
2635
2636 log("read_pixels", "Reading data . . .\n");
2637
2638 i=-1;
2639 icount = 0;
2640
2641 while ( ! qefile.eof() ) {
2642
2643 // get line from the file
2644
2645 qefile.getline(line, LINE_MAX_LENGTH);
2646
2647 // skip if comment
2648
2649 if ( *line == '#' )
2650 continue;
2651
2652 // if it is the first valid value, it is the number of QE data points
2653
2654 if ( i < 0 ) {
2655
2656 // get the number of datapoints
2657
2658 sscanf(line, "%d", &pointsQE);
2659
2660 // allocate memory for the table of QEs
2661
2662 QE = new float ** [ct_NPixels];
2663
2664 for ( i=0; i<ct_NPixels; ++i ) {
2665 QE[i] = new float * [2];
2666 QE[i][0] = new float[pointsQE];
2667 QE[i][1] = new float[pointsQE];
2668 }
2669
2670 QElambda = new float [pointsQE];
2671
2672 for ( i=0; i<pointsQE; ++i ) {
2673 qefile.getline(line, LINE_MAX_LENGTH);
2674 sscanf(line, "%f", &QElambda[i]);
2675 }
2676
2677 i=0;
2678
2679 continue;
2680 }
2681
2682 // get the values (num-pixel, num-datapoint, QE-value)
2683
2684 if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2685 break;
2686
2687 if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2688 ((j-1) < pointsQE) && ((j-1) > -1) ) {
2689 QE[i-1][0][j-1] = QElambda[j-1];
2690 QE[i-1][1][j-1] = qe;
2691 }
2692
2693 if ( i > ct_NPixels) break;
2694
2695 icount++;
2696
2697 }
2698
2699 if(icount/pointsQE < ct_NPixels){
2700 error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2701 icount/pointsQE, ct_NPixels );
2702 }
2703
2704 // close file
2705
2706 qefile.close();
2707
2708 // test QE
2709
2710 for(icount=0; icount< ct_NPixels; icount++){
2711 for(i=0; i<pointsQE; i++){
2712 if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
2713 QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
2714 error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
2715 icount, i, QE[icount][0][i], QE[icount][1][i] );
2716 }
2717 }
2718 }
2719
2720 // end
2721
2722 log("read_pixels", "Done.\n");
2723
2724}
2725//!@}
2726
2727
2728//!-----------------------------------------------------------
2729// @name pixels_are_neig
2730//
2731// @desc check whether two pixels are neighbours
2732//
2733// @var pix1 Number of the first pixel
2734// @var pix2 Number of the second pixel
2735// @return TRUE: both pixels are neighbours; FALSE: oth.
2736//
2737// @date Wed Sep 9 17:58:37 MET DST 1998
2738//------------------------------------------------------------
2739// @function
2740
2741//!@{
2742int
2743pixels_are_neig(int pix1, int pix2)
2744{
2745 if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2746 SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2747 > ct_PixelWidth_corner_2_corner )
2748 return ( FALSE );
2749 else
2750 return ( TRUE );
2751}
2752//!@}
2753
2754//!-----------------------------------------------------------
2755// @name igen_pixel_coordinates
2756//
2757// @desc generate the pixel center coordinates
2758//
2759// @var *pcam structure camera containing all the
2760// camera information
2761// @return total number of pixels
2762//
2763// DP
2764//
2765// @date Thu Oct 14 10:41:03 CEST 1999
2766//------------------------------------------------------------
2767// @function
2768
2769//!@{
2770/******** igen_pixel_coordinates() *********************************/
2771
2772int igen_pixel_coordinates(struct camera *pcam) {
2773 /* generate pixel coordinates, return value is number of pixels */
2774
2775 int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2776 float fsegment_fract;
2777 double dtsize;
2778 double dhsize;
2779 double dpsize;
2780 double dxfirst_pix;
2781 double dyfirst_pix;
2782 double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5;
2783 double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5;
2784
2785
2786 double dstartx, dstarty; /* for the gap pixels and outer pixels */
2787 int j, nrow;
2788
2789 dpsize = pcam->dpixdiameter_cm;
2790 dtsize = dpsize * sqrt(3.) / 2.;
2791 dhsize = dpsize / 2.;
2792
2793 /* Loop over central pixels to generate co-ordinates */
2794
2795 for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2796
2797 /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2798
2799 pcam->dpixsizefactor[ipixno-1] = 1.;
2800
2801 in = 0;
2802
2803 i = 0;
2804 itot_inside_ring = 0;
2805 iring_no = 0;
2806
2807 /* Calculate the number of pixels out to and including the ring containing pixel number */
2808 /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2809
2810 while (itot_inside_ring == 0){
2811
2812 iN = 3*(i*(i+1)) + 1;
2813
2814 if (ipixno <= iN){
2815 iring_no = i;
2816 itot_inside_ring = iN;
2817 }
2818
2819 i++;
2820 }
2821
2822
2823 /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2824
2825 ipix_in_ring = 0;
2826 for (i = 0; i < iring_no; ++i){
2827
2828 ipix_in_ring = ipix_in_ring + 6;
2829 }
2830
2831 /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2832 /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2833 /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2834
2835 isegment = 0;
2836 fsegment_fract = 0.;
2837 if (iring_no > 0) {
2838
2839 isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2840
2841 fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2842
2843 }
2844
2845 /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2846 /* pixel width (flat to flat) dpsize. */
2847
2848 dxfirst_pix = dpsize*iring_no;
2849 dyfirst_pix = 0.;
2850
2851 /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2852
2853 ddxseg1 = - dhsize*iring_no;
2854 ddyseg1 = dtsize*iring_no;
2855 ddxseg2 = -dpsize*iring_no;
2856 ddyseg2 = 0.;
2857 ddxseg3 = ddxseg1;
2858 ddyseg3 = -ddyseg1;
2859 ddxseg4 = -ddxseg1;
2860 ddyseg4 = -ddyseg1;
2861 ddxseg5 = -ddxseg2;
2862 ddyseg5 = 0.;
2863
2864 /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2865 /* anti-clockwise around the ring by adding the segment to segment vectors. */
2866
2867 switch (isegment) {
2868
2869 case 0:
2870
2871 pcam->dxc[ipixno-1] = 0.;
2872 pcam->dyc[ipixno-1] = 0.;
2873
2874 case 1:
2875 pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2876 pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2877
2878 break;
2879
2880 case 2:
2881
2882 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2883 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2884
2885 break;
2886
2887 case 3:
2888
2889 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
2890 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
2891
2892 break;
2893
2894 case 4:
2895
2896 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
2897 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
2898
2899 break;
2900
2901 case 5:
2902
2903 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
2904 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2905
2906 break;
2907
2908 case 6:
2909
2910 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2911 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2912
2913 break;
2914
2915 default:
2916
2917 fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2918 return(0);
2919
2920 } /* end switch */
2921
2922 } /* end for */
2923
2924 dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2925 dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2926
2927 if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2928
2929 j = pcam->inumcentralpixels;
2930
2931 for(i=0; i<pcam->inumgappixels; i=i+6){
2932 pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2933 pcam->dyc[j + i ] = dstarty;
2934 pcam->dpixsizefactor[j + i] = 1.;
2935 pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2936 pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2937 pcam->dpixsizefactor[j + i + 1] = 1.;
2938 pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2939 pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2940 pcam->dpixsizefactor[j + i+ 2] = 1.;
2941 pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2942 pcam->dyc[j + i + 3] = dstarty;
2943 pcam->dpixsizefactor[j + i+ 3] = 1.;
2944 pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2945 pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2946 pcam->dpixsizefactor[j + i+ 4] = 1.;
2947 pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2948 pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2949 pcam->dpixsizefactor[j + i + 5] = 1.;
2950 } /* end for */
2951 } /* end if */
2952
2953 /* generate positions of the outer pixels */
2954
2955 if( pcam->inumbigpixels > 0 ){
2956
2957 j = pcam->inumcentralpixels + pcam->inumgappixels;
2958
2959 for(i=0; i<pcam->inumbigpixels; i++){
2960 pcam->dpixsizefactor[j + i] = 2.;
2961 }
2962
2963 in = 0;
2964
2965 nrow = (int) ceil(dstartx / 2. / dpsize);
2966
2967 while(in < pcam->inumbigpixels){
2968
2969 pcam->dxc[j + in] = dstartx + dpsize;
2970 pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2971 pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2972 pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2973 pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2974 pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2975 pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2976 pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2977 pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2978 pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2979 pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2980 pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2981 for(i=1; i<nrow; i++){
2982 pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2983 pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2984 pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2985 pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2986 pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2987 pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2988 pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2989 pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2990 pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
2991 pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
2992 pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
2993 pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
2994 }
2995 in = in + 6 * nrow;
2996 dstartx = dstartx + 2. * dpsize;
2997 nrow = nrow + 1;
2998 } /* end while */
2999
3000 } /* end if */
3001
3002 return(pcam->inumpixels);
3003
3004}
3005//!@}
3006
3007//!-----------------------------------------------------------
3008// @name bpoint_is_in_pix
3009//
3010// @desc check if a point (x,y) in camera coordinates is inside a given pixel
3011//
3012// @var *pcam structure camera containing all the
3013// camera information
3014// @var dx, dy point coordinates in centimeters
3015// @var ipixnum pixel number (starting at 0)
3016// @return TRUE if the point is inside the pixel, FALSE otherwise
3017//
3018// DP
3019//
3020// @date Thu Oct 14 16:59:04 CEST 1999
3021//------------------------------------------------------------
3022// @function
3023
3024//!@{
3025
3026/******** bpoint_is_in_pix() ***************************************/
3027
3028#define sqrt13 0.577350269 // = 1./sqrt(3.)
3029#define sqrt3 1.732050807 // = sqrt(3.)
3030
3031int bpoint_is_in_pix(double dx, double dy, MGeomCamMagic *pgeomcam)
3032{
3033 /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
3034 /* the pixel is assumed to be a "closed set" */
3035
3036 /*
3037 a = length of one of the edges of one pixel,
3038 b = half the width of one pixel
3039 */
3040
3041 const int numN = pgeomcam->GetNumPixels();
3042
3043 for (int i=0; i<numN; i++)
3044 {
3045 MGeomPix &pixel = (*pgeomcam)[i];
3046 const double b = pixel.GetR()/2;
3047 const double a = pixel.GetR()/sqrt3;
3048
3049 const double xx = dx - pixel.GetX();
3050 const double yy = dy - pixel.GetY();
3051
3052 if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
3053 ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){
3054
3055 return i; // inside i
3056 }
3057
3058 // return -1; // outside
3059 }
3060
3061 return -1; // outside
3062}
3063
3064//!@}
3065
3066//------------------------------------------------------------
3067// @name dist_r_P
3068//
3069// @desc distance straight line r - point P
3070//
3071// @date Sat Jun 27 05:58:56 MET DST 1998
3072// @function @code
3073//------------------------------------------------------------
3074// dist_r_P
3075//
3076// distance straight line r - point P
3077//------------------------------------------------------------
3078
3079float
3080dist_r_P(float a, float b, float c,
3081 float l, float m, float n,
3082 float x, float y, float z)
3083{
3084 return (
3085 sqrt((SQR((a-x)*m-(b-y)*l) +
3086 SQR((b-y)*n-(c-z)*m) +
3087 SQR((c-z)*l-(a-x)*n))/
3088 (SQR(l)+SQR(m)+SQR(n))
3089 )
3090 );
3091}
3092
3093//------------------------------------------------------------
3094// @name check_reflector_file
3095//
3096// @desc check if a given reflector file has the right signature
3097// @desc return TRUE or FALSE
3098//
3099// @date Mon Feb 14 16:44:21 CET 2000
3100// @function @code
3101//------------------------------------------------------------
3102
3103int check_reflector_file(FILE *infile){
3104
3105 char sign[20]; // auxiliary variable
3106
3107 strcpy(sign, REFL_SIGNATURE_B);
3108
3109 fread( (char *)sign, strlen(REFL_SIGNATURE_B), 1, infile);
3110 if (strcmp(sign, REFL_SIGNATURE_A) == 0){
3111 fread( (char *)sign, 1, 1, infile);
3112 return 4;
3113 }
3114 else if (strcmp(sign, REFL_SIGNATURE_B) == 0){
3115 fread( (char *)sign, 1, 1, infile);
3116 return 5;
3117 }
3118 else if (strcmp(sign, REFL_SIGNATURE_C) == 0){
3119 fread( (char *)sign, 1, 1, infile);
3120 return 6;
3121 }
3122 else {
3123 cout << "ERROR: Signature of .rfl file is not correct\n";
3124 cout << '"' << sign << '"' << '\n';
3125 cout << "should be: " << REFL_SIGNATURE_A <<" or "<< REFL_SIGNATURE_B <<" or " << REFL_SIGNATURE_C <<" or "<< '\n';
3126 return(FALSE);
3127 }
3128
3129
3130}
3131
3132//------------------------------------------------------------
3133// @name lin_interpol
3134//
3135// @desc interpolate linearly between two points returning the
3136// @desc y-value of the result
3137//
3138// @date Thu Feb 17 11:31:32 CET 2000
3139// @function @code
3140//------------------------------------------------------------
3141
3142float lin_interpol( float x1, float y1, float x2, float y2, float x){
3143
3144 if( (x2 - x1)<=0. ){ // avoid division by zero, return average
3145 cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
3146 return((y1+y2)/2.);
3147 }
3148 else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
3149 return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
3150 }
3151}
3152
3153
3154//------------------------------------------------------------
3155// @name produce_phes
3156//
3157// @desc read the photons of an event, pixelize them and simulate QE
3158// @desc return various statistics and the array of Photoelectrons
3159//
3160// @date Mon Feb 14 16:44:21 CET 2000
3161// @function @code
3162//------------------------------------------------------------
3163
3164int produce_phes( FILE *sp, // the input file
3165 MGeomCamMagic *camgeom, // the camera layout
3166 float minwl_nm, // the minimum accepted wavelength
3167 float maxwl_nm, // the maximum accepted wavelength
3168 class MTrigger *trigger, // the generated phes
3169 class MFadc *fadc,
3170 int *itotnphe, // total number of produced photoelectrons
3171 float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
3172 int *incph, // total number of cph read
3173 float *tmin_ns, // minimum arrival time of all phes
3174 float *tmax_ns, // maximum arrival time of all phes
3175 int star // star=0 starfield rotation star=-1 not rotation
3176 ){
3177
3178 static uint i;
3179 static int k, ipixnum;
3180 static float cx, cy, wl, qe, t;
3181 static float cu, cv, cw;
3182 static MCCphoton photon;
3183 static float **qept;
3184 static char flag[SIZE_OF_FLAGS + 1];
3185 static float radius_mm;
3186 static float SFR , C3 , C2 , C1;
3187 const double pi = 3.14159;
3188
3189 // reset variables
3190
3191 for ( i=0; i<camgeom->GetNumPixels(); ++i ){
3192
3193 nphe[i] = 0.0;
3194
3195 }
3196
3197 *incph = 0;
3198
3199 radius_mm = camgeom->GetMaxRadius();
3200
3201
3202 //- - - - - - - - - - - - - - - - - - - - - - - - -
3203 // read photons and "map" them into the pixels
3204 //--------------------------------------------------
3205
3206 // initialize CPhoton
3207
3208 photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
3209
3210 // read the photons data
3211
3212
3213 // loop over the photons
3214
3215 while (1) {
3216
3217 fread ( flag, SIZE_OF_FLAGS, 1, sp );
3218
3219 if (isA( flag, FLAG_END_OF_EVENT ))
3220 break;
3221
3222 memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
3223
3224 fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
3225
3226
3227 // increase number of photons
3228
3229 (*incph)++;
3230
3231 // Chceck if photon is inside trigger time range
3232
3233 t = photon.get_t() ;
3234
3235 if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
3236 continue;
3237 }
3238
3239 //
3240 // Pixelization
3241 //
3242
3243 // Where the photon provide(0 from starfield, -1 from Mmcs ???
3244
3245 if ( star == 0 && Starfield_rotate == TRUE )
3246
3247 {
3248 get_teles_axis(&zenith, &azimutal);
3249 // Introduction SFR angle
3250
3251 zenith = (pi/180) * zenith;
3252
3253 azimutal =( pi/180) * azimutal;
3254
3255 C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
3256
3257 C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
3258
3259 C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
3260
3261 SFR = acos( C1/C2 );
3262
3263 if ( sin(azimutal) < 0)
3264 {
3265 SFR = 2 * pi - SFR;
3266 }
3267 else
3268 {
3269 SFR = SFR;
3270 }
3271
3272//final of introduction SFR angle//
3273
3274 cx = photon.get_x() * cos(SFR) - photon.get_y() * sin(SFR);
3275 cy = photon.get_y() * cos(SFR) + photon.get_x() * sin(SFR);
3276
3277 }
3278 else
3279 {
3280 cx = photon.get_x();
3281 cy = photon.get_y();
3282
3283 }
3284
3285
3286 // get wavelength
3287
3288 wl = photon.get_wl();
3289
3290 // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n";
3291
3292 // check if photon has valid wavelength and is inside outermost camera radius
3293
3294 if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy)*10 > radius_mm ) ){
3295 continue;
3296
3297 }
3298
3299 ipixnum = bpoint_is_in_pix(cx*10, cy*10, camgeom);
3300
3301 // -1 = the photon is in none of the pixels
3302 // 0 = the phton is in the central pixel, which is not used for trigger
3303 if (ipixnum==-1 || ipixnum==0) {
3304 continue;
3305 }
3306
3307 //+++
3308 // QE simulation
3309 //---
3310
3311 // set pointer to the QE table of the relevant pixel
3312
3313 qept = (float **)QE[ipixnum];
3314
3315 // check if wl is inside table; outside the table, QE is assumed to be zero
3316
3317 if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
3318 continue;
3319
3320 }
3321
3322 // find data point in the QE table (-> k)
3323
3324 k = 1; // start at 1 because the condition was already tested for 0
3325 while (k < pointsQE-1 && qept[0][k] < wl){
3326 k++;
3327 }
3328
3329 // calculate the qe value between 0. and 1.
3330
3331 qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
3332
3333 // Apply incient angular correction due to Light Guides
3334
3335 cu=photon.get_u();
3336 cv=photon.get_v();
3337 cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*pi/180.0;
3338
3339 // find data point in the QE table (-> k)
3340
3341 k = 0;
3342 while (k < pointsWC-1 && WC[0][k] < cw){
3343 k++;
3344 }
3345
3346 // correct the qe with WC data.
3347
3348 qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
3349
3350 // if random > quantum efficiency, reject it
3351
3352 if ( (RandomNumber) > qe ) {
3353 continue;
3354 }
3355
3356 //+++
3357 // The photon has produced a photo electron
3358 //---
3359
3360 // cout << " accepted\n";
3361
3362 // increment the number of photoelectrons in the relevant pixel
3363
3364 nphe[ipixnum] += 1.0;
3365
3366 // store the new photoelectron
3367
3368 fadc->Fill(ipixnum,(t-*tmin_ns) , trigger->FillShow(ipixnum,t-*tmin_ns));
3369
3370 *itotnphe += 1;
3371 }
3372
3373 return(0);
3374
3375}
3376
3377
3378//------------------------------------------------------------
3379// @name produce_nsbrates
3380//
3381// @desc read the starfield file, call produce_phes on it in,
3382// @desc different wavebands, calculate the nsbrates
3383//
3384// @date Mon Feb 14 16:44:21 CET 2000
3385// @function @code
3386//------------------------------------------------------------
3387
3388int produce_nsbrates( char *iname, // the starfield input file name
3389 MGeomCamMagic *camgeom, // camera layout
3390 float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
3391 // the NSB rates in phe/ns for each pixel
3392 ){
3393
3394 uint i, j; // counters
3395 int k, ii; // counters
3396
3397 MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
3398 MFadc flashadc;
3399
3400 static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
3401 WAVEBANDBOUND2,
3402 WAVEBANDBOUND3,
3403 WAVEBANDBOUND4,
3404 WAVEBANDBOUND5,
3405 WAVEBANDBOUND6 };
3406
3407 FILE *infile; // the input file
3408 fpos_t fileposition; // marker on the input file
3409 static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
3410 static MCEventHeader evth; // the event header
3411 static MCEventHeader evth_2; // the event header
3412 static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
3413 int reflector_file_version;
3414 int itnphe; // total number of produced photoelectrons
3415 int itotnphe; // total number of produced photoelectrons after averaging
3416 int incph; // total number of cph read
3417 float tmin_ns; // minimum arrival time of all phes
3418 float tmax_ns; // maximum arrival time of all phes
3419 float integtime_ns; // integration time from the starfield generator
3420
3421 // open input file
3422
3423 log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
3424
3425 infile = fopen( iname, "r" );
3426
3427 // check if the satrfield input file exists
3428
3429 if ( infile == NULL ) {
3430
3431 log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
3432
3433 log( SIGNATURE, "There is not NSB from the Stars\n");
3434
3435 return (0);
3436 }
3437
3438 // get signature, and check it
3439
3440 if((reflector_file_version=check_reflector_file( infile ))==FALSE){
3441 exit(1);
3442 }
3443
3444 // initialize flag
3445
3446 strcpy( cflag, " \0" );
3447
3448 // get flag
3449
3450 fread( cflag, SIZE_OF_FLAGS, 1, infile );
3451
3452 if ( ! feof(infile)){
3453
3454 // reading .rfl file
3455
3456 if(!isA( cflag, FLAG_START_OF_RUN )){
3457 error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
3458 }
3459 else { // found start of run
3460
3461 fread( cflag, SIZE_OF_FLAGS, 1, infile );
3462
3463 if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
3464
3465 // get MCEventHeader
3466
3467 if (reflector_file_version<6)
3468 fread( (char*)&evth, evth.mysize(), 1, infile );
3469 else
3470 fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
3471
3472 if (reflector_file_version<6)
3473 integtime_ns = evth.get_energy();
3474 else
3475 integtime_ns = evth_2.get_energy();
3476
3477 // memorize where we are in the file
3478
3479 if (fgetpos( infile, &fileposition ) != 0){
3480 error( SIGNATURE, "Cannot position in file ...\n");
3481 }
3482
3483 // loop over the wavebands
3484
3485 for(i=0; i<iNUMWAVEBANDS; i++){
3486
3487 // initialize the rate array
3488
3489 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
3490 rate_phepns[j][i] = 0.;
3491 }
3492
3493 itotnphe = 0;
3494
3495 // read the photons and produce the photoelectrons
3496 // - in order to average over the QE simulation, call the
3497 // production function iNUMNSBPRODCALLS times
3498
3499 for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
3500
3501 // position the file pointer to the beginning of the photons
3502
3503 fsetpos( infile, &fileposition);
3504
3505 // produce photoelectrons
3506
3507 k = produce_phes( infile,
3508 camgeom,
3509 wl_nm[i],
3510 wl_nm[i+1],
3511 &trigger, // this is a dummy here
3512 &flashadc, // this is a dummy here
3513 &itnphe,
3514 nphe, // we want this!
3515 &incph,
3516 &tmin_ns,
3517 &tmax_ns,
3518 0 ); // photons from starfield
3519
3520 if( k != 0 ){ // non-zero returnvalue means error
3521 cout << "Exiting.\n";
3522 exit(1);
3523 }
3524
3525 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
3526 rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
3527 }
3528
3529 itotnphe += itnphe;
3530
3531 } // end for(ii=0 ...
3532
3533 fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
3534 wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
3535
3536 } // end for(i=0 ...
3537
3538 }
3539 else{
3540 cout << "Starfield file contains no event.\nExiting.\n";
3541 exit(1);
3542 } // end if( isA ... event
3543 } // end if ( isA ... run
3544 }
3545 else{
3546 cout << "Starfield file contains no run.\nExiting.\n";
3547 exit(1);
3548 }
3549
3550 fclose( infile );
3551
3552 return(0);
3553
3554}
3555
3556
3557//------------------------------------------------------------
3558// @name produce_nsb_phes
3559//
3560// @desc produce the photoelectrons from the nsbrates
3561//
3562// @date Thu Feb 17 17:10:40 CET 2000
3563// @function @code
3564//------------------------------------------------------------
3565
3566int produce_nsb_phes( float *atmin_ns,
3567 float *atmax_ns,
3568 float theta_rad,
3569 struct camera *cam,
3570 float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
3571 float difnsbr_phepns[iMAXNUMPIX],
3572 float extinction[iNUMWAVEBANDS],
3573 float fnpx[iMAXNUMPIX],
3574 MTrigger *trigger,
3575 MFadc *fadc,
3576 int *inphe,
3577 float base_mv[iMAXNUMPIX]){
3578
3579 float simtime_ns; // NSB simulation time
3580 int i, j, k, ii;
3581 float zenfactor; // correction factor calculated from the extinction
3582 int inumnsbphe; // number of photoelectrons caused by NSB
3583
3584 float t;
3585 TRandom random; // Random numbers generator
3586
3587 random.SetSeed ((UInt_t) (RandomNumber*1000));
3588
3589 ii = *inphe; // avoid dereferencing
3590
3591 // check if the arrival times are set; if not generate them
3592
3593 if(*atmin_ns <SIMTIMEOFFSET_NS || *atmin_ns > *atmax_ns){
3594 *atmin_ns = 0.;
3595 *atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
3596
3597 }
3598 else{ // extend the event time window by the given offsets
3599
3600 *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
3601 *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
3602
3603 simtime_ns = *atmax_ns - *atmin_ns;
3604
3605 // make sure the simulated time is long enough for the FADC
3606 // simulation and not too long
3607
3608 if(simtime_ns< TOTAL_TRIGGER_TIME){
3609 *atmin_ns = *atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
3610 *atmax_ns = *atmin_ns + TOTAL_TRIGGER_TIME;
3611 simtime_ns = TOTAL_TRIGGER_TIME;
3612 }
3613
3614 if(simtime_ns> TOTAL_TRIGGER_TIME){
3615 *atmax_ns =*atmin_ns + TOTAL_TRIGGER_TIME;
3616 simtime_ns = TOTAL_TRIGGER_TIME;
3617 }
3618
3619 }
3620
3621 // initialize baselines
3622
3623 for(i=0; i<cam->inumpixels; i++){
3624 base_mv[i] = 0.;
3625 }
3626
3627 // calculate baselines and generate phes
3628
3629 for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
3630
3631 // calculate the effect of the atmospheric extinction
3632
3633 zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
3634
3635 for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
3636
3637 inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
3638 * simtime_ns );
3639
3640 base_mv[j] += inumnsbphe;
3641
3642 // randomize
3643
3644 if (inumnsbphe>0.0){
3645 inumnsbphe = random.Poisson (inumnsbphe );
3646 }
3647
3648 // create the photoelectrons
3649
3650 for(k=0; k < inumnsbphe; k++){
3651
3652 t=(RandomNumber * simtime_ns);
3653
3654 (*fadc).Fill(j,t ,(*trigger).FillNSB(j,t));
3655
3656 ii++; // increment total number of photoelectons
3657
3658 fnpx[j] += 1.; // increment number of photoelectrons in this pixel
3659
3660 }
3661
3662 } // end for(j=0 ...
3663 } // end for(i=0 ...
3664
3665 // finish baseline calculation
3666
3667 for(i=0; i<cam->inumpixels; i++){
3668 base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
3669 }
3670
3671 *inphe = ii; // update the pointer
3672
3673 return(0);
3674}
3675
3676
3677// @endcode
3678
3679
3680//=------------------------------------------------------------
3681//!@subsection Log of this file.
3682
3683//!@{
3684//
3685// $Log: not supported by cvs2svn $
3686// Revision 1.45 2002/10/29 17:15:27 blanch
3687// Reading several reflector versions.
3688//
3689// Revision 1.44 2002/10/18 16:53:03 blanch
3690// Modification to read several reflector files.
3691//
3692// Revision 1.43 2002/09/13 10:53:39 blanch
3693// Minor change to remove some undisired comments.
3694//
3695// Revision 1.42 2002/09/09 16:00:49 blanch
3696// Statement has been included to avoid writting to disk MParContainer and MArray.
3697// It has also been added the effect of the WC, the actual values must be added,
3698// once they are measured.
3699//
3700// Revision 1.41 2002/09/04 09:57:42 blanch
3701// Modifications done to use MGeomCam from MARS.
3702//
3703// Revision 1.40 2002/07/16 16:15:22 blanch
3704// A first implementation of the Star field rotation has been done.
3705//
3706// Revision 1.39 2002/04/27 10:48:39 blanch
3707// Some unused varibles have been removed.
3708//
3709// Revision 1.38 2002/03/18 18:44:29 blanch
3710// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
3711//
3712// Revision 1.37 2002/03/18 16:42:20 blanch
3713// The data member fProductionSite of the MMcRunHeader has been set to 0,
3714// which means that the prodution site is unknown.
3715//
3716// Revision 1.35 2002/03/15 15:15:52 blanch
3717// Several modifications to simulate the actual trigger zone.
3718//
3719// Revision 1.34 2002/03/13 18:13:56 blanch
3720// Some changes to fill correctly the new format of MMcRunHeader.
3721//
3722// Revision 1.33 2002/03/04 17:21:48 blanch
3723// Small and not important changes.
3724//
3725// Revision 1.32 2002/02/28 15:04:52 blanch
3726// A small back has been solved. Before, while not using the option
3727// writte_all_images, not all triggered showers were stored. Now it is solved.
3728// For that it is importatn taht the less restrictive trigger option is
3729// checked first.
3730// A new facility has been introduced and now one can choose the step size in
3731// trigger loop mode for the discriminator threshold.
3732// The close-compact topology for less than 3 pixels does not make sense. Before
3733// the program was ignoring that, now it switch to simple neighbour condition.
3734//
3735// Revision 1.31 2002/01/18 17:41:02 blanch
3736// The option of adding noise to all pixels or to not adding the noise
3737// has been added.
3738// We removed the pixels larger than 577. When there were more than one
3739// trigger in one shower, the pixel number was increasing. Now it is
3740// flagged by the variable MMcTrig::fFirstLvlTrig.
3741//
3742// Revision 1.30 2001/11/27 09:49:54 blanch
3743// Fixing bug which was treating wrongly the extension of star photons.
3744//
3745// Revision 1.29 2001/11/14 17:38:23 blanch
3746// Sveral changes have been done:
3747// - bpoint_in_in_pixel has been dodified to speed up the program
3748// - Some minor changes have been done to remove warnings
3749// - buffer size and split version of the Branches have been removed
3750// - Some modifications were needed to adat the program to the new
3751// MRawEvtData::DeletePixels
3752//
3753// Revision 1.28 2001/10/26 16:31:45 blanch
3754// Removing several warnings.
3755//
3756// Revision 1.27 2001/09/05 10:04:33 blanch
3757// *** empty log message ***
3758//
3759// Revision 1.26 2001/07/19 09:29:53 blanch
3760// Different threshold for each pixel can be used.
3761//
3762// Revision 1.25 2001/05/08 08:07:54 blanch
3763// New numbering for branches from different trigger conditions has been
3764// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
3765// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
3766// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
3767// information already filled.
3768//
3769// Revision 1.24 2001/04/06 16:48:09 magicsol
3770// New camera version able to read the new format of the reflector output:
3771// reflector 0.4
3772//
3773// Revision 1.23 2001/03/28 16:13:41 blanch
3774// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
3775// been solved.
3776//
3777// Revision 1.22 2001/03/20 18:52:43 blanch
3778// *** empty log message ***
3779//
3780// Revision 1.21 2001/03/19 19:53:03 blanch
3781// Some print outs have been removed.
3782//
3783// Revision 1.20 2001/03/19 19:30:06 magicsol
3784// Minor changes have been done to improve the FADC pedestals treatment.
3785//
3786// Revision 1.19 2001/03/05 11:14:41 magicsol
3787// I changed the position of readinf a parameter. It is a minnor change.
3788//
3789// Revision 1.18 2001/03/05 10:36:52 blanch
3790// A branch with information about the FADC simulation (MMcFadcHeader) is writen
3791// in the McHeaders tree of the aoutput root file.
3792// The NSB contribution is only applied if the the number of phe form the shower
3793// are above a given number.
3794//
3795// Revision 1.17 2001/02/23 11:05:57 magicsol
3796// Small changes due to slightly different output format and the introduction of
3797// pedesals for teh FADC.
3798//
3799// Revision 1.16 2001/01/18 18:44:40 magicsol
3800// *** empty log message ***
3801//
3802// Revision 1.15 2001/01/17 09:32:27 harald
3803// The changes are neccessary to have the same name for trees in MC and in
3804// data. So now there should be no differences in MC data and real data from
3805// FADC system.
3806//
3807// Revision 1.14 2001/01/15 12:33:34 magicsol
3808// Some modifications have been done to use the new (Dec'2000) Raw data format.
3809// There are also some minnors modifications to adapt the improvements in the
3810// MTrigger class (overlaping time and trigger cells).
3811//
3812// Revision 1.13 2000/10/25 08:14:23 magicsol
3813// The routine that produce poisson random numbers to decide how many phe
3814// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
3815// is used.
3816//
3817// Revision 1.12 2000/09/22 17:40:18 harald
3818// Added a lot of changes done by oscar.
3819//
3820// Revision 1.11 2000/09/21 11:47:33 harald
3821// Oscar found some smaller errors in the calculation of the pixel shape and
3822// corrected it.
3823//
3824// $Log: not supported by cvs2svn $
3825// Revision 1.45 2002/10/29 17:15:27 blanch
3826// Reading several reflector versions.
3827//
3828// Revision 1.44 2002/10/18 16:53:03 blanch
3829// Modification to read several reflector files.
3830//
3831// Revision 1.43 2002/09/13 10:53:39 blanch
3832// Minor change to remove some undisired comments.
3833//
3834// Revision 1.42 2002/09/09 16:00:49 blanch
3835// Statement has been included to avoid writting to disk MParContainer and MArray.
3836// It has also been added the effect of the WC, the actual values must be added,
3837// once they are measured.
3838//
3839// Revision 1.41 2002/09/04 09:57:42 blanch
3840// Modifications done to use MGeomCam from MARS.
3841//
3842// Revision 1.40 2002/07/16 16:15:22 blanch
3843// A first implementation of the Star field rotation has been done.
3844//
3845// Revision 1.39 2002/04/27 10:48:39 blanch
3846// Some unused varibles have been removed.
3847//
3848// Revision 1.38 2002/03/18 18:44:29 blanch
3849// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
3850//
3851// Revision 1.37 2002/03/18 16:42:20 blanch
3852// The data member fProductionSite of the MMcRunHeader has been set to 0,
3853// which means that the prodution site is unknown.
3854//
3855// Revision 1.35 2002/03/15 15:15:52 blanch
3856// Several modifications to simulate the actual trigger zone.
3857//
3858// Revision 1.34 2002/03/13 18:13:56 blanch
3859// Some changes to fill correctly the new format of MMcRunHeader.
3860//
3861// Revision 1.33 2002/03/04 17:21:48 blanch
3862// Small and not important changes.
3863//
3864// Revision 1.32 2002/02/28 15:04:52 blanch
3865// A small back has been solved. Before, while not using the option
3866// writte_all_images, not all triggered showers were stored. Now it is solved.
3867// For that it is importatn taht the less restrictive trigger option is
3868// checked first.
3869// A new facility has been introduced and now one can choose the step size in
3870// trigger loop mode for the discriminator threshold.
3871// The close-compact topology for less than 3 pixels does not make sense. Before
3872// the program was ignoring that, now it switch to simple neighbour condition.
3873//
3874// Revision 1.31 2002/01/18 17:41:02 blanch
3875// The option of adding noise to all pixels or to not adding the noise
3876// has been added.
3877// We removed the pixels larger than 577. When there were more than one
3878// trigger in one shower, the pixel number was increasing. Now it is
3879// flagged by the variable MMcTrig::fFirstLvlTrig.
3880//
3881// Revision 1.30 2001/11/27 09:49:54 blanch
3882// Fixing bug which was treating wrongly the extension of star photons.
3883//
3884// Revision 1.29 2001/11/14 17:38:23 blanch
3885// Sveral changes have been done:
3886// - bpoint_in_in_pixel has been dodified to speed up the program
3887// - Some minor changes have been done to remove warnings
3888// - buffer size and split version of the Branches have been removed
3889// - Some modifications were needed to adat the program to the new
3890// MRawEvtData::DeletePixels
3891//
3892// Revision 1.28 2001/10/26 16:31:45 blanch
3893// Removing several warnings.
3894//
3895// Revision 1.27 2001/09/05 10:04:33 blanch
3896// *** empty log message ***
3897//
3898// Revision 1.26 2001/07/19 09:29:53 blanch
3899// Different threshold for each pixel can be used.
3900//
3901// Revision 1.25 2001/05/08 08:07:54 blanch
3902// New numbering for branches from different trigger conditions has been
3903// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
3904// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
3905// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
3906// information already filled.
3907//
3908// Revision 1.24 2001/04/06 16:48:09 magicsol
3909// New camera version able to read the new format of the reflector output:
3910// reflector 0.4
3911//
3912// Revision 1.23 2001/03/28 16:13:41 blanch
3913// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
3914// been solved.
3915//
3916// Revision 1.22 2001/03/20 18:52:43 blanch
3917// *** empty log message ***
3918//
3919// Revision 1.21 2001/03/19 19:53:03 blanch
3920// Some print outs have been removed.
3921//
3922// Revision 1.20 2001/03/19 19:30:06 magicsol
3923// Minor changes have been done to improve the FADC pedestals treatment.
3924//
3925// Revision 1.19 2001/03/05 11:14:41 magicsol
3926// I changed the position of readinf a parameter. It is a minnor change.
3927//
3928// Revision 1.18 2001/03/05 10:36:52 blanch
3929// A branch with information about the FADC simulation (MMcFadcHeader) is writen
3930// in the McHeaders tree of the aoutput root file.
3931// The NSB contribution is only applied if the the number of phe form the shower
3932// are above a given number.
3933//
3934// Revision 1.17 2001/02/23 11:05:57 magicsol
3935// Small changes due to slightly different output format and the introduction of
3936// pedesals for teh FADC.
3937//
3938// Revision 1.16 2001/01/18 18:44:40 magicsol
3939// *** empty log message ***
3940//
3941// Revision 1.15 2001/01/17 09:32:27 harald
3942// The changes are neccessary to have the same name for trees in MC and in
3943// data. So now there should be no differences in MC data and real data from
3944// FADC system.
3945//
3946// Revision 1.14 2001/01/15 12:33:34 magicsol
3947// Some modifications have been done to use the new (Dec'2000) Raw data format.
3948// There are also some minnors modifications to adapt the improvements in the
3949// MTrigger class (overlaping time and trigger cells).
3950//
3951// Revision 1.13 2000/10/25 08:14:23 magicsol
3952// The routine that produce poisson random numbers to decide how many phe
3953// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
3954// is used.
3955//
3956// Revision 1.10 2000/07/04 14:10:20 MagicSol
3957// Some changes have been done in the root output file. The RawEvt tree is only
3958// stored in the single trigger mode.
3959// The trigger input parameters are also given by the general input card.
3960// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
3961// each one can be studied seperately.
3962//
3963// Revision 1.9 2000/06/13 13:25:24 blanch
3964// The multiple arrays have been replaced, since they do not work
3965// in alpha machines. Now we are using pointers and the command new
3966// to allocate memory.
3967//
3968// Revision 1.8 2000/05/11 13:57:27 blanch
3969// The option to loop over several trigger configurations has been included.
3970// Some non-sense with triggertime range has been solved.
3971// Montecarlo information and ADC counts are saved in a root file.
3972// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
3973//
3974// Revision 1.7 2000/03/24 18:10:46 blanch
3975// A first FADC simulation and a trigger simulation are already implemented.
3976// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
3977// 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.
3978//
3979// Revision 1.6 2000/03/20 18:35:11 blanch
3980// 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.
3981//
3982// Revision 1.5 2000/02/18 17:40:35 petry
3983// This version includes drastic changes compared to camera.cxx 1.4.
3984// It is not yet finished and not immediately useful because the
3985// trigger simulation is not yet re-implemented. I had to take it
3986// out together with some other stuff in order to tidy the whole
3987// program up. This is not meant as an insult to anyone. I needed
3988// to do this in order to be able to work on it.
3989//
3990// This version has been put in the repository in order to be
3991// able to share the further development with others.
3992//
3993// If you need something working, wait or take an earlier one.
3994// See file README.
3995//
3996// Revision 1.4 2000/01/25 08:36:23 petry
3997// The pixelization in previous versions was buggy.
3998// This is the first version with a correct pixelization.
3999//
4000// Revision 1.3 2000/01/20 18:22:17 petry
4001// Found little bug which makes camera crash if it finds a photon
4002// of invalid wavelength. This bug is now fixed and the range
4003// of valid wavelengths extended to 290 - 800 nm.
4004// This is in preparation for the NSB simulation to come.
4005// Dirk
4006//
4007// Revision 1.2 1999/11/19 08:40:42 harald
4008// Now it is possible to compile the camera programm under osf1.
4009//
4010// Revision 1.1.1.1 1999/11/05 11:59:31 harald
4011// This the starting point for CVS controlled further developments of the
4012// camera program. The program was originally written by Jose Carlos.
4013// But here you can find a "rootified" version to the program. This means
4014// that there is no hbook stuff in it now. Also the output of the
4015// program changed to the MagicRawDataFormat.
4016//
4017// The "rootification" was done by Dirk Petry and Harald Kornmayer.
4018//
4019// Revision 1.3 1999/10/22 15:01:28 petry
4020// version sent to H.K. and N.M. on Fri Oct 22 1999
4021//
4022// Revision 1.2 1999/10/22 09:44:23 petry
4023// first synthesized version which compiles and runs without crashing;
4024//
4025// Revision 1.1.1.1 1999/10/21 16:35:10 petry
4026// first synthesised version
4027//
4028// Revision 1.13 1999/03/15 14:59:05 gonzalez
4029// camera-1_1
4030//
4031// Revision 1.12 1999/03/02 09:56:10 gonzalez
4032// *** empty log message ***
4033//
4034//
4035//!@}
4036
4037//=EOF
Note: See TracBrowser for help on using the repository browser.