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

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