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

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