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

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