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

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