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

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