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

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