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

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