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

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