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

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