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

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