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

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