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

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