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

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