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

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