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

Last change on this file since 4308 was 3044, checked in by blanch, 21 years ago
Compiling with root 3.05 and updating MARS files.
File size: 129.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.68 $
24// $Author: blanch $
25// $Date: 2004-02-06 17:39:24 $
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#include "TObjArray.h"
61
62#include "MTrigger.hxx"
63#include "MFadc.hxx"
64#include "MLons.hxx"
65
66#include "MRawRunHeader.h"
67#include "MRawEvtData.h"
68#include "MRawEvtHeader.h"
69#include "MRawCrateArray.h"
70#include "MRawCrateData.h"
71#include "MMcEvt.hxx"
72#include "MMcTrig.hxx"
73#include "MMcRunHeader.hxx"
74#include "MMcConfigRunHeader.h"
75#include "MMcCorsikaRunHeader.h"
76#include "MMcTrigHeader.hxx"
77#include "MMcFadcHeader.hxx"
78#include "MGeomCamMagic.h"
79#include "MGeomCamMagic919.h"
80#include "MGeomCamMagicHG.h"
81#include "MGeomCamECO1000.h"
82#include "MGeomCamECO1000HG.h"
83#include "MGeomCamCT1.h"
84#include "MGeomCamCT1Daniel.h"
85#include "MGeomPix.h"
86
87/*!@"
88
89 All the defines are located in the file |camera.h|.
90
91 @"*/
92
93#include "camera.h"
94
95//!@}
96
97/*!@"
98
99 The following set of flags are used in time of compilation. They do
100 not affect directly the behaviour of the program at run-time
101 (though, of course, if you disconnected the option for
102 implementation of the Trigger logic, you will not be able to use any
103 trigger at all. The 'default' values mean default in the sense of
104 what you got from the server when you obtained this program.
105
106 @"*/
107
108//!@{
109
110// flag for debugging (default: OFF )
111#define __DEBUG__
112#undef __DEBUG__
113
114
115//!@}
116
117//=-----------------------------------------------------------
118//!@subsection Definition of global variables.
119
120/*!@"
121
122 Now we define some global variables with data about the telescope,
123 such as "focal distance", number of pixels/mirrors,
124 "size of the camera", and so on.
125
126 @"*/
127
128/*!@"
129
130 Depending on the telescope we are using (CT1 or MAGIC), the
131 information stored in the definition file is different.
132 The variable |ct_Type| has the value 0 when we use
133 CT1, and 1 when we use MAGIC.
134
135 @"*/
136
137/*!@"
138
139 And this is the information about the whole telescope.
140
141 @"*/
142
143//!@{
144
145// parameters of the CT (from the CT definition file)
146
147//@: Number of pixels
148static int ct_NPixels;
149
150//@: Number of CT
151static int ct_Number;
152
153//@: Camera geometries
154static int ct_Geometry;
155
156//@: list of showers to be skipped
157static int *Skip;
158
159//@: number of showers to be skipped
160static int nSkip=0;
161
162//@: flag: TRUE: data come from STDIN; FALSE: from file
163static int Data_From_STDIN = FALSE;
164
165//@: flag: TRUE: write all images to output; FALSE: only triggered showers
166static int Write_All_Images = FALSE;
167
168static int Write_McEvt = TRUE;
169static int Write_McTrig = TRUE;
170static int Write_McFADC = TRUE;
171static int Write_RawEvt = FALSE;
172
173//@: flag: TRUE: selection on the energy
174static int Select_Energy = TRUE;
175
176//@: Lower edge of the selected energy range (in GeV)
177static float Select_Energy_le = 0.0;
178
179//@: Upper edge of the selected energy range (in GeV)
180static float Select_Energy_ue = 100000.0;
181
182//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
183static int FADC_Scan = FALSE;
184
185//@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't
186static int Trigger_Scan = FALSE;
187
188//@: flag: TRUE: loop trigger analysis over several thresholds, multiplicities and topologies; FALSE: a single trigger configuration
189static int Trigger_Loop = FALSE;
190
191//@: flag: TRUE: Different threshold for each pixel ; FALSE: same threshold for all pixels
192static int Individual_Thres_Pixel = FALSE;
193
194//@: Properties of the trigger
195static float Trigger_gate_length = 6.0;
196static float Trigger_response_ampl = 1.0;
197static float Trigger_response_fwhm = 2.0;
198static float Trigger_overlaping_time= 0.25;
199static float Trigger_noise= 0.3;
200
201//@: Properties of the FADC
202static float FADC_response_ampl = MFADC_RESPONSE_INTEGRAL;
203static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
204static float FADC_resp_ampl_out = MFADC_RESPONSE_INTEGRAL;
205static float FADC_resp_fwhm_out = MFADC_RESPONSE_FWHM;
206static float FADC_noise = 2.0;
207static float DIGITAL_noise = 2.0;
208static float FADC_high2low = HIGH2LOWGAIN;
209
210//@: Trigger conditions for a single trigger mode
211static float **qThreshold;
212static int Trigger_multiplicity[100];
213static int Trigger_topology[100];
214
215//@: Upper and lower edges of the trigger loop
216static float Trigger_loop_lthres = 2.0;
217static float Trigger_loop_uthres = 10.0;
218static float Trigger_loop_sthres = 1.0;
219static int Trigger_loop_lmult = 2;
220static int Trigger_loop_umult = 10;
221static int Trigger_loop_ltop = 0;
222static int Trigger_loop_utop = 2;
223
224//!@}
225
226/*!@"
227
228 The following double-pointer is a 2-dimensional table with information
229 about each pixel. The routine read_pixels will generate
230 the information for filling it using igen_pixel_coordinates().
231
232 @"*/
233
234//!@{
235// Pointer to a tables/Arrays with information about the pixels
236// and data stored on them with information about the pixels
237
238
239//@: contents of the pixels (ph.e.)
240static float *fnpix;
241
242
243//!@}
244
245/*!@"
246
247 The following double-pointer is a 2-dimensional table with the
248 Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
249 of the wavelength @$\lambda@$. The routine |read_pixels()| will read
250 also this information from the file |qe.dat|.
251
252 @"*/
253
254//!@{
255// Pointer to a table with QE, number of datapoints, and wavelengths
256
257//@: table of QE
258static float ****QE;
259
260//@: number of datapoints for the QE curve
261static int pointsQE[100];
262
263//@: table of QE
264static float *QElambda;
265
266//@: table of lightguide = WC;
267static float **WC;
268
269//@: number of datapoints for the WC curve
270static int pointsWC;
271
272//!@}
273
274/*!@"
275
276 The following double-pointer is a 2-dimensional table with information
277 about each mirror in the dish. The routine |read_ct_file()| will read
278 this information from the CT definition file.
279
280 @"*/
281
282/*!@"
283
284 We define a table into where random numbers will be stored.
285 The routines used for random number generation are provided by
286 |RANLIB| (taken from NETLIB, |www.netlib.org|), and by
287 the routine |double drand48(void)| (prototype defined in
288 |stdlib.h|) through the macro |RandomNumber| defined in
289 |camera.h|.
290
291 @"*/
292
293
294//!@}
295
296extern char FileName[];
297// switch on starfield rotation
298static int Starfield_rotate = TRUE;
299
300//=-----------------------------------------------------------
301// @subsection Main program.
302
303//!@{
304
305//++++++++++++++++++++++++++++++++++++++++
306// MAIN PROGRAM
307//----------------------------------------
308
309int main(int argc, char **argv)
310{
311
312 //!@' @#### Definition of variables.
313 //@'
314
315 char inname_CT[100][256]; //@< array of names for each CT input file
316 char starfieldname[256]; //@< starfield input file name
317 char qe_filename[256]; //@< qe input file name
318
319 char datname[256]; //@< data (ASCII) output file name
320
321 char rootname[256] ; //@< ROOT file name
322 char rootname_loop[256] ; //@< ROOT file name
323
324 char parname[256]; //@< parameters file name
325
326 char nsbpathname[256]; //@< directory with the dataabse for th NSB
327 char nsbpath_outer[256]; //@< directory with the dataabse for outer pixels
328
329 char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
330 char flag_new[SIZE_OF_FLAGS + 1]; //@< New flag for the run header in the .rfl file
331 char GeometryName[100][50]; //@< Name of MGeomCam for each CT
332 int GeometryCamera[100]; //@< Identification of MGeomCam for each CT
333 int TriggerPixels[100]; //@< Number of pixels in the trigger region for each CT
334
335 int reflector_file_version=0; //@< vrsion of he reflector file
336
337 FILE *inputfile[100]; //@< stream for the input file
338
339 ofstream datafile; //@< stream for the data file
340
341 MCRunHeader mcrunh; //@< Run Header class (MC)
342 MCEventHeader mcevth[100]; //@< Event Header class (MC)
343 MCEventHeader_2 mcevth_2[100]; //@< Event Header class (MC) for reflector > 0.6
344 MCCphoton cphoton; //@< Cherenkov Photon class (MC)
345
346 int inumphe; //@< number of photoelectrons in an event from showers
347 int inumphe_CT[100]; //@< number of photoelectrons in an event from showers
348 float inumphensb[100]; //@< number of photoelectrons in an event from nsb
349
350 float arrtmin_ns; //@ arrival time of the first photoelectron
351 float arrtmax_ns; //@ arrival time of the last photoelectron
352
353 float thetaCT[100], phiCT[100]; //@< parameters of a given shower
354 float thetashw, phishw; //@< parameters of a given shower
355 float coreX, coreY; //@< core position
356 float impactD[100]; //@< impact parameter
357 float l1, m1, n1; //@< auxiliary variables
358 float l2, m2, n2; //@< auxiliary variables
359 float num, den; //@< auxiliary variables
360 float factorqe_NSB[100]; //@< factor on the NSB depending of QE
361
362 int nshow=0; //@< partial number of shower in a given run
363 int ntshow=0; //@< total number of showers
364 int ncph[100]; //@< partial number of photons in a given run
365 int ntcph[100]; //@< total number of photons
366
367 int ibr, j, k; //@< simple counters
368
369 int addElecNoise; //@< Will we add ElecNoise?
370
371 float riseDiskThres;
372 float secureDiskThres;
373
374 int simulateNSB; //@< Will we simulate NSB?
375 int nphe2NSB; //@< From how many phe will we simulate NSB?
376 float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel)
377 float diffnsb_phepns[100][iMAXNUMPIX]; //@< diffuse NSB values for each pixel
378 //@< derived from meanNSB
379 float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS]; //@< non-diffuse nsb
380 //@< photoelectron rates
381 float nsb_phepns[100][iMAXNUMPIX]; //@< NSB values for each pixel
382 float nsb_phepns_rotated[100][iMAXNUMPIX]; //@< NSB values for each pixel after rotation
383
384 Float_t nsb_trigresp[TRIGGER_TIME_SLICES]; //@< array to write the trigger
385 //@< response from the database
386 Float_t nsb_fadcresp[(Int_t) SLICES_MFADC]; //@< array to write the fadc
387 //@< response from the database
388 Byte_t trigger_map[((Int_t)(CAMERA_PIXELS/8))+1]; //@< Pixels on when the
389 //@< camera triggers
390 Float_t fadc_elecnoise[CAMERA_PIXELS]; //@< Electronic noise for each pixel
391 Float_t fadc_diginoise[CAMERA_PIXELS]; //@< Digital noise for each pixel
392 Float_t fadc_pedestals[CAMERA_PIXELS]; //@< array for fadc pedestals values
393 Float_t fadc_sigma[CAMERA_PIXELS]; //@< array for fadc pedestals sigma
394 Float_t fadc_sigma_low[CAMERA_PIXELS]; //@< array for fadc pedestals sigma
395
396 float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
397 EXTWAVEBAND1,
398 EXTWAVEBAND2,
399 EXTWAVEBAND3,
400 EXTWAVEBAND4,
401 EXTWAVEBAND5
402 };
403 float zenfactor=1.0; // correction factor calculated from the extinction
404
405 int nstoredevents = 0;
406 int flagstoring = 0;
407
408 int ntrigger[100]; //@< number of triggers in the whole file
409 int btrigger = 0; //@< trigger flag
410 int ithrescount; //@< counter for loop over threshold trigger
411 float fthrescount; //@< value for loop over threshold trigger
412 int imulticount; //@< counter for loop over multiplicity trigger
413 int itopocount; //@< counter for loop over topology trigger
414 int isorttopo[3]; //@< sorting the topologies
415 int icontrigger; //@< number of trigger conditions to be analised
416
417 float fpixelthres[CAMERA_PIXELS]; //@< Threshold values
418
419 TArrayC *fadcValues; //@< the analog Fadc High gain signal for pixels
420 TArrayC *fadcValuesLow; //@< the analog Fadc Low gain signal for pixels
421
422 int still_in_loop = FALSE;
423
424 //@< Variables to fill the McRunHeader
425 Int_t sfRaH = 5;
426 Int_t sfRaM = 34;
427 Int_t sfRaS = 32;
428 Int_t sfDeD = 22;
429 Int_t sfDeM = 00;
430 Int_t sfDeS = 52;
431 Float_t shthetamax = 0.0;
432 Float_t shthetamin = 0.0;
433 Float_t shphimax = 0.0;
434 Float_t shphimin = 0.0;
435 Float_t telestheta = 0.0;
436 Float_t telesphi = 0.0;
437 Float_t sofftheta = 0.0;
438 Float_t soffphi = 0.0;
439 UInt_t corsika = 5200 ;
440 Float_t maxpimpact = 0.0;
441
442 // Star Field Rotation variables
443 Float_t zenith = 0.0;
444 Float_t azimutal = 90.0;
445 Float_t rho , C3 , C2 , C1;
446
447 //!@' @#### Definition of variables for |getopt()|.
448 //@'
449
450 int ch, errflg = 0; //@< used by getopt
451
452 //!@' @#### Initialisation of some variables
453 //@'
454
455 qThreshold = new float *[100];
456 for(int i=0;i<100;i++)
457 qThreshold[i] = new float [CAMERA_PIXELS];
458
459 for(int i=0;i<iMAXNUMPIX;i++){
460 for(int ict=0;ict<100;ict++){
461 nsb_phepns[ict][i]=0;
462 nsb_phepns_rotated[ict][i]=0;
463 }
464 for(j=0;j<iNUMWAVEBANDS;j++)
465 nsbrate_phepns[i][j]=0.0; //@< Starlight rates initialised at 0.0
466 }
467 for(int i=0;i<100;i++)
468 ntcph[i]=0;
469 /*!@'
470
471 @#### Beginning of the program.
472
473 We start with the main program. First we (could) make some
474 presentation, and follows the reading of the parameters file (now
475 from the |stdin|), the reading of the CT parameters file, and the
476 creation of the output file, where the processed data will be
477 stored.
478
479 */
480
481 //++
482 // START
483 //--
484
485 // make unbuffered output
486
487 // cout.setf ( ios::stdio );
488
489 // parse command line options (see reflector.h)
490
491 parname[0] = '\0';
492
493 optarg = NULL;
494 while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
495 switch (ch) {
496 case 'f':
497 strcpy(parname, optarg);
498 break;
499 case 'h':
500 usage();
501 break;
502 default :
503 errflg++;
504 }
505
506 // show help if error
507
508 if ( errflg>0 )
509 usage();
510
511 // make some sort of presentation
512
513 present();
514
515 // read parameters file
516
517 if ( strlen(parname) < 1 )
518 readparam(NULL);
519 else
520 readparam(parname);
521
522 // read data from file or from STDIN?
523
524 Data_From_STDIN = get_data_from_stdin();
525
526 // get number of telescopes and camera geometries
527
528 ct_Number=get_ct_number();
529 ct_Geometry=get_ct_geometry();
530
531 for(int ict=0;ict<ct_Number;ict++){
532 ntrigger[ict]=0;
533 if(ct_Geometry>=ict*10){
534 GeometryCamera[ict]=int(ct_Geometry/pow(10.0,ict))%10;
535 }
536 else
537 GeometryCamera[ict]=0;
538 }
539
540 // structure holding the camera definition
541
542 TObjArray camgeom;
543
544 for (int i=0; i<ct_Number;i++){
545 switch(GeometryCamera[i]){
546 case 1: camgeom[i]=new MGeomCamMagic;
547 strcpy(GeometryName[i],"MGeomCamMagic");
548 TriggerPixels[i]=TRIGGER_PIXELS_1;
549 break;
550 case 2: camgeom[i]=new MGeomCamMagic919;
551 strcpy(GeometryName[i],"MGeomCamMagic919");
552 TriggerPixels[i]=TRIGGER_PIXELS_2;
553 break;
554 case 3: camgeom[i]=new MGeomCamMagicHG;
555 strcpy(GeometryName[i],"MGeomCamMagicHG");
556 TriggerPixels[i]=TRIGGER_PIXELS_3;
557 break;
558 case 5: camgeom[i]=new MGeomCamECO1000;
559 strcpy(GeometryName[i],"MGeomCamECO1000");
560 TriggerPixels[i]=TRIGGER_PIXELS_5;
561 break;
562 case 6: camgeom[i]=new MGeomCamECO1000HG;
563 strcpy(GeometryName[i],"MGeomCamECO1000HG");
564 TriggerPixels[i]=TRIGGER_PIXELS_6;
565 break;
566 case 8: camgeom[i]=new MGeomCamCT1;
567 strcpy(GeometryName[i],"MGeomCamCT1");
568 TriggerPixels[i]=TRIGGER_PIXELS_8;
569 break;
570 case 9: camgeom[i]=new MGeomCamCT1Daniel;
571 strcpy(GeometryName[i],"MGeomCamCT1Daniel");
572 TriggerPixels[i]=TRIGGER_PIXELS_9;
573 break;
574 default:
575 error( SIGNATURE, "Camera geometry %i is not defined. Stoping camera program ...", GeometryCamera[i]);
576 exit(1);
577 break;
578 }
579 }
580
581 ct_NPixels=0;
582
583 for(int ict=0;ict<ct_Number;ict++)
584 ct_NPixels=(((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels()>(UInt_t)ct_NPixels)?
585 (Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels():
586 ct_NPixels;
587
588 // read WC and QE files
589
590 // FIX ME !
591 // Currently the PMT N of any camera with same average QE has the same QE.
592
593 QE = new float *** [ct_Number];
594 for(int ict=0;ict<ct_Number;ict++){
595 strcpy( qe_filename, get_qe_filename(ict));
596 read_QE(qe_filename,ict);
597 }
598
599 read_WC();
600
601 // write all images, even those without trigger?
602
603 Write_All_Images = get_write_all_events();
604
605 Write_McEvt = get_write_McEvt() ;
606 Write_McTrig = get_write_McTrig() ;
607 Write_McFADC = get_write_McFadc() ;
608 Write_RawEvt = get_write_RawEvt() ;
609
610 FADC_Scan = get_FADC_Scan();
611 Trigger_Scan = get_Trigger_Scan();
612
613 Individual_Thres_Pixel = get_indi_thres_pixel();
614
615 get_secure_threhold(&riseDiskThres,&secureDiskThres);
616
617 get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm,
618 &FADC_resp_ampl_out, &FADC_resp_fwhm_out);
619 FADC_high2low=get_High_to_Low();
620
621 // FIXME --- tirgger properties may depend on the Camrea geometry
622
623 get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);
624
625 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);
626
627 icontrigger =((int)((Trigger_loop_uthres-Trigger_loop_lthres)
628 /Trigger_loop_sthres)+1)*
629 (Trigger_loop_umult-Trigger_loop_lmult+1)*
630 (Trigger_loop_utop-Trigger_loop_ltop+1);
631
632 // Trigger loop operation is not implemented for Multi telscopes
633
634 if ( Trigger_Loop && ct_Number > 1 ){
635 cout<<"ERROR : camera::main : Trigger loop option is not";
636 cout<<" implemented for Multi Telescopes option. "<<icontrigger<<
637 " "<<ct_Number<<endl;
638 exit(1);
639 }
640
641 if (!Trigger_Loop){
642 get_Trigger_Single (qThreshold,
643 &Trigger_multiplicity[0],
644 &Trigger_topology[0]);
645 icontrigger=1;
646 }
647 else
648 get_threshold(qThreshold[0]);
649
650 // get filenames
651
652 for(int ict=0;ict<ct_Number;ict++)
653 strcpy( inname_CT[ict], get_input_filename(ict) );
654
655 strcpy( starfieldname, get_starfield_filename() );
656 strcpy( datname, get_data_filename() );
657 strcpy( rootname, get_root_filename() );
658 strcpy( rootname_loop, get_loop_filename() );
659 strcpy( nsbpathname, get_nsb_directory() );
660 strcpy( nsbpath_outer, get_nsb_directory_outer() );
661
662 // get different parameters of the simulation
663
664 addElecNoise = add_elec_noise(&FADC_noise, &DIGITAL_noise, &Trigger_noise);
665 simulateNSB = get_nsb( &meanNSB, &nphe2NSB );
666
667
668 // get selections on the parameters
669
670 Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
671
672 //Parameters for starfield rotation
673
674 Starfield_rotate = get_starfield_rotate();
675
676 // log filenames information
677 for(Int_t ict=0;ict<ct_Number;ict++){
678
679 log(SIGNATURE,"\t%s : %i\n\t%20s:\t%s\n",
680 "------- TELESCOPE ",ict+1,
681 "Geometry ", GeometryName[ict]);
682 strcpy( qe_filename, get_qe_filename(ict));
683
684 // Look to factor for NSB respect to emi_coat PMTs
685 if(strstr(qe_filename,"qe-emi-coat.RFL.dat") != 0)
686 factorqe_NSB[ict]=EMICOAT_NSB;
687 else if(strstr(qe_filename,"qe-emi.RFL.dat") != 0)
688 factorqe_NSB[ict]=EMI_NSB;
689 else if(strstr(qe_filename, "qe-intevac_hpd.RFL.dat") != 0)
690 factorqe_NSB[ict]=HPD_NSB;
691 else{
692 log(SIGNATURE,"WARNING : the factor for the diffuse NSB for this QE file (%s) is not known. The number of phe(ns stated in the inputcard will be used.\n",qe_filename);
693 factorqe_NSB[ict]=1.0;
694 }
695 log(SIGNATURE,
696 "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s\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",
697 "Filenames",
698 "In", inname_CT[ict],
699 "Stars", starfieldname,
700 "NSB database","(inner pixels)", nsbpathname,
701 "(outer pixels)", nsbpath_outer,
702 "QE", qe_filename,
703 "Data", datname,
704 "ROOT", rootname
705 );
706
707 // log parameters information
708
709 log(SIGNATURE,
710 "%s:\n\t%20s: %f: %s\n\t%20s: %f\n\t%20s: %f\n",
711 "Parameters",
712 "NSB (phes/ ns 0.1*0.1 deg^2 239 m^2)", meanNSB*factorqe_NSB[ict], ONoff(simulateNSB),
713 "Pedestals = ", get_FADC_pedestal(),
714 "High to Low gain = ", FADC_high2low
715 );
716
717 // log Trigger information
718
719 if (Trigger_Loop) {
720 log(SIGNATURE,
721 "%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",
722 "Trigger Loop mode",
723 "Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres,
724 "Multiplicity",Trigger_loop_lmult,Trigger_loop_umult,
725 "Topology",Trigger_loop_ltop,Trigger_loop_utop,
726 rootname_loop);
727 }
728 else if (Individual_Thres_Pixel == FALSE){
729 log(SIGNATURE,
730 "%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n",
731 "Single Trigger mode",
732 "Threshold",qThreshold[ict][0],
733 "Multiplicity",Trigger_multiplicity[ict],
734 "Topology",Trigger_topology[ict]);
735 }
736 else{
737 log(SIGNATURE,
738 "%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n",
739 "Single Trigger mode",
740 "Threshold","Different for each pixel",
741 "Multiplicity",Trigger_multiplicity[ict],
742 "Topology",Trigger_topology[ict]);
743 }
744 log(SIGNATURE,"\t%s\n",
745 "END TELESCOPE --------");
746 }
747 // log flags information
748
749 log(SIGNATURE,
750 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %3.2f(%s) %3.2f+%3.2f(%s) %s\n",
751 "Flags",
752 "Data_From_STDIN", ONoff(Data_From_STDIN),
753 "Write_All_Events", ONoff(Write_All_Images),
754 "Rotate Starfield", ONoff(Starfield_rotate),
755 "Electronic Noise", Trigger_noise,"trigger",
756 FADC_noise,DIGITAL_noise,"fadc",ONoff(addElecNoise)
757 );
758
759 // log flags information
760
761 log(SIGNATURE,
762 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
763 "Root output",
764 "Write_McEvt", ONoff(Write_McEvt),
765 "Write_McTrig", ONoff(Write_McTrig),
766 "Write_McFADC", ONoff(Write_McFADC),
767 "Write_RawEvt", ONoff(Write_RawEvt));
768
769 // log selections
770
771 log(SIGNATURE,
772 "%s:\n\t%20s: %s (%f:%f)\n",
773 "Selections:",
774 "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
775
776 // Definition and initialization of array to save trigger statistics
777
778 int ***ntriggerloop;
779
780 ntriggerloop= new int ** [(int)((Trigger_loop_uthres-Trigger_loop_lthres)
781 /Trigger_loop_sthres)];
782 for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;fthrescount+=Trigger_loop_sthres, ithrescount++){
783 ntriggerloop[ithrescount]= new int * [Trigger_loop_umult-Trigger_loop_lmult+1];
784 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
785 ntriggerloop[ithrescount][imulticount]= new int [Trigger_loop_utop-Trigger_loop_ltop+1];
786 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
787 ntriggerloop[ithrescount][imulticount][itopocount]=0;
788 }
789 }
790 }
791
792 // We should be careful that topologies are sort from
793 // the less to the more restrictive one.
794
795 if (Trigger_loop_utop==Trigger_loop_ltop)
796 for(int is=0; is<3;is++)
797 isorttopo[is]=is;
798 else {
799 isorttopo[0]=1;
800 isorttopo[1]=0;
801 isorttopo[2]=2;
802 }
803
804 // get list of showers to evt. skip
805
806 nSkip = get_nskip_showers();
807
808 if (nSkip > 0) {
809 Skip = new int[ nSkip ];
810 get_skip_showers( Skip );
811
812 log(SIGNATURE, "There are some showers to skip:\n");
813 for (int i=0; i<nSkip; ++i)
814 log(SIGNATURE, "\tshower # %d\n", Skip[i]);
815 }
816
817
818 // read parameters from the ct.def file
819
820 //read_ct_file();
821
822 Int_t Lev0, Lev1;
823 Int_t Lev0MT[100], Lev1MT[100];
824
825 fadcValues = new TArrayC(FADC_SLICES);
826 fadcValuesLow = new TArrayC(FADC_SLICES);
827
828 // number of pixels for parameters
829
830 // Switched off writing TObject
831
832 MArray::Class()->IgnoreTObjectStreamer();
833 MParContainer::Class()->IgnoreTObjectStreamer();
834
835 // initialise ROOT
836
837 TROOT simple("simple", "MAGIC Telescope Monte Carlo");
838
839 // initialise instance of Trigger and FADC classes
840
841 MTrigger **Trigger_CT;
842 Trigger_CT = new MTrigger *[ct_Number];
843
844 for (int i=0; i<ct_Number;i++){
845 Trigger_CT[i] = new MTrigger(TriggerPixels[i],
846 ((MGeomCam*)(camgeom.UncheckedAt(i))),
847 Trigger_gate_length,
848 Trigger_overlaping_time,
849 Trigger_response_ampl,
850 Trigger_response_fwhm); //@< A instance of the Class MTrigger
851 Trigger_CT[i]->SetSeed(UInt_t(i+get_seeds(0)));
852 }
853
854 for (int ict=0; ict<ct_Number;ict++){
855 // Generate database for trigger electronic noise
856 if (addElecNoise)
857 Trigger_CT[ict]->SetElecNoise(Trigger_noise) ;
858
859 // Set Right Discriminator threshold, taking into account trigger pixels
860 Trigger_CT[ict]->CheckThreshold(&qThreshold[ict][0],GeometryCamera[ict]);
861 }
862 // Set flag in pixel 0 (not used for trigger) that indicates if secure pixel
863 // is active: secureDiskThres*10000+riseDiskThres
864
865 for (int ict=0; ict<ct_Number;ict++){
866 if(riseDiskThres>0.0)
867 qThreshold[ict][0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres;
868 else
869 qThreshold[ict][0]=0.0;
870 }
871 // Initialise McTrig information class if we want to save trigger informtion
872
873 MMcTrig **McTrig = NULL;
874 MMcTrigHeader **HeaderTrig = NULL;
875 MMcFadcHeader **HeaderFadc = NULL;
876
877 int numberBranches;
878
879 if (!Trigger_Loop)
880 numberBranches=ct_Number;
881 else
882 numberBranches=icontrigger;
883
884 if (Write_McTrig){
885
886 McTrig = new MMcTrig * [numberBranches];
887
888 for (int i=0;i<numberBranches;i++) {
889 McTrig[i] = new MMcTrig();
890 }
891
892 HeaderTrig = new MMcTrigHeader * [numberBranches];
893
894 for (int i=0;i<numberBranches;i++) {
895 HeaderTrig[i] = new MMcTrigHeader();
896 }
897 }
898
899 if (Write_McFADC){
900
901 HeaderFadc = new MMcFadcHeader * [numberBranches];
902
903 for (int i=0;i<numberBranches;i++) {
904 HeaderFadc[i] = new MMcFadcHeader();
905 }
906 }
907
908 MFadc **Fadc_CT;
909 Fadc_CT = new MFadc *[ct_Number];
910
911 for (int i=0; i<ct_Number;i++)
912 Fadc_CT[i] =
913 new MFadc(((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels(),
914 FADC_response_ampl,FADC_response_fwhm,
915 FADC_resp_ampl_out,FADC_resp_fwhm_out,
916 get_trig_delay()
917 ) ; //@< A instance of the Class MFadc
918
919
920 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
921 //
922 // Set the FADC pedestals for that run
923 // Some modifications
924 // mut be done to simulate a more realistic distribution of the pedestals.
925 // This simualtion is done int the SetPedestals methode inside the
926 // class MFadc
927 // Currentlly a input_pedestal array is declared with the pedestals.
928 // Thy can also be set randomly following a flat distribution.
929 //
930 /////////////////////////////////////////////////////////////////////
931
932 Float_t input_pedestals[ct_NPixels];
933
934 for(int i=0;i<ct_NPixels;i++)
935 input_pedestals[i]=get_FADC_pedestal();
936 for (int ict=0; ict<ct_Number;ict++){
937
938 Fadc_CT[ict]->SetPedestals(input_pedestals);
939 Fadc_CT[ict]->SetHigh2LowGain(FADC_high2low);
940
941 // Generate database for the Fadc electronic noise
942 if (addElecNoise){
943 Fadc_CT[ict]->SetElecNoise(FADC_noise);
944 Fadc_CT[ict]->SetDigitalNoise(DIGITAL_noise);
945 }
946 }
947
948 // Prepare the raw data output
949
950 // Header Tree
951 MGeomCam **camdummy = new MGeomCam * [ct_Number];
952 for (int ict=0; ict<ct_Number;ict++){
953 camdummy[ict]= ((MGeomCam*)(camgeom.UncheckedAt(ict)));
954 }
955 MRawRunHeader *RunHeader= new MRawRunHeader();
956 MMcRunHeader *McRunHeader = new MMcRunHeader();
957 MMcCorsikaRunHeader *McCorsikaRunHeader = new MMcCorsikaRunHeader();
958
959 MMcConfigRunHeader **McConfigRunHeader = NULL;
960 McConfigRunHeader = new MMcConfigRunHeader * [numberBranches];
961 for (int i=0;i<numberBranches;i++) {
962 McConfigRunHeader[i] = new MMcConfigRunHeader();
963 }
964
965 // Header branch
966
967 MRawEvtHeader **EvtHeader = NULL;
968
969 if (Write_RawEvt) {
970 EvtHeader = new MRawEvtHeader * [numberBranches];
971
972 for (int i=0;i<numberBranches;i++) {
973 EvtHeader[i] = new MRawEvtHeader();
974 }
975 }
976
977 // Data branch
978
979 MRawEvtData **EvtData = NULL; // Data branch
980
981 if (Write_RawEvt) {
982 EvtData = new MRawEvtData * [numberBranches];
983
984 for (int i=0;i<numberBranches;i++) {
985 EvtData[i] = new MRawEvtData();
986 EvtData[i]->Init(RunHeader); // We need thr RunHeader to read
987 // number of pixels
988 }
989 }
990
991 MMcEvt **McEvt = NULL;
992
993 if (Write_McEvt) {
994 if (ct_Number>1){
995 McEvt = new MMcEvt *[numberBranches];
996 for (int i=0;i<numberBranches;i++)
997 McEvt[i]=new MMcEvt();
998 }
999 else {
1000 McEvt = new MMcEvt *[1];
1001 McEvt[0]=new MMcEvt();
1002 }
1003 }
1004 //
1005 // initalize a temporal ROOT file
1006 //
1007
1008 TFile outfile_temp ( rootname , "RECREATE" );
1009
1010 // create a Tree for the Header Event
1011 TTree HeaderTree("RunHeaders","Headers of Run");
1012
1013 // define branches of Header Tree
1014
1015 char help[4];
1016
1017 HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
1018 &RunHeader);
1019
1020 HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
1021 &McRunHeader);
1022
1023 HeaderTree.Branch("MMcCorsikaRunHeader","MMcCorsikaRunHeader",
1024 &McCorsikaRunHeader);
1025
1026 if(!Trigger_Loop && Write_McTrig && ct_Number==1){
1027
1028 HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader",
1029 &HeaderTrig[0]);
1030 }
1031 if (ct_Number==1){
1032 HeaderTree.Branch("MGeomCam",GeometryName[0],
1033 &camdummy[0]);
1034 HeaderTree.Branch("MMcConfigRunHeader","MMcConfigRunHeader",
1035 &McConfigRunHeader[0]);
1036 }
1037 else{
1038 char branchname[20];
1039 for (int ict=0; ict<ct_Number;ict++){
1040 sprintf(help,"%i",ict+1);
1041 strcpy (branchname, "MGeomCam;");
1042 strcat (branchname, & help[0]);
1043 strcat (branchname, ".");
1044 HeaderTree.Branch(branchname,GeometryName[ict],
1045 &camdummy[ict]);
1046 }
1047 for(int i=0;i<numberBranches;i++){
1048 sprintf(help,"%i",i+1);
1049 strcpy (branchname, "MMcConfigRunHeader;");
1050 strcat (branchname, & help[0]);
1051 strcat (branchname, ".");
1052 HeaderTree.Branch(branchname,"MMcConfigRunHeader",
1053 &McConfigRunHeader[i]);
1054 }
1055
1056 }
1057
1058 if ((Trigger_Loop || ct_Number>1) && Write_McTrig){
1059 ibr=0;
1060 for(char branchname[20];ibr<numberBranches;ibr++){
1061
1062 sprintf(help,"%i",ibr+1);
1063 strcpy (branchname, "MMcTrigHeader;");
1064 strcat (branchname, & help[0]);
1065 strcat (branchname, ".");
1066 HeaderTree.Branch(branchname,"MMcTrigHeader",
1067 &HeaderTrig[ibr]);
1068 }
1069 }
1070
1071 if(!Trigger_Loop && Write_McFADC && ct_Number==1){
1072
1073 HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader",
1074 &HeaderFadc[0]);
1075 }
1076 if ((Trigger_Loop || ct_Number>1) && Write_McFADC){
1077 ibr=0;
1078 for(char branchname[20];ibr<numberBranches;ibr++){
1079
1080 sprintf(help,"%i",ibr+1);
1081 strcpy (branchname, "MMcFadcHeader;");
1082 strcat (branchname, & help[0]);
1083 strcat (branchname, ".");
1084 HeaderTree.Branch(branchname,"MMcFadcHeader",
1085 &HeaderFadc[ibr]);
1086 }
1087 }
1088
1089 // Fill branches for MRawRunHeader
1090
1091 RunHeader->SetMagicNumber(kMagicNumber);
1092 RunHeader->SetFormatVersion(4);
1093 RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
1094 RunHeader->SetRunType(256);
1095 RunHeader->SetRunNumber(0);
1096 RunHeader->SetNumSamples(0,FADC_SLICES);
1097 RunHeader->SetNumCrates(1);
1098 RunHeader->SetNumPixInCrate(ct_NPixels);
1099
1100 // Fill branches for MMcTrigHeader
1101
1102 if(!Trigger_Loop && Write_McTrig && ct_Number==1){
1103
1104 HeaderTrig[0]->SetTopology((Short_t) Trigger_topology[0]);
1105 HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity[0]);
1106 HeaderTrig[0]->SetThreshold(qThreshold[0]);
1107 HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
1108 HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
1109 HeaderTrig[0]->SetOverlap(Trigger_overlaping_time);
1110 HeaderTrig[0]->SetGate(Trigger_gate_length);
1111 HeaderTrig[0]->SetElecNoise(Trigger_noise);
1112 }
1113 if(!Trigger_Loop && Write_McTrig && ct_Number>1){
1114 for(int i=0;i<ct_Number;i++){
1115 HeaderTrig[i]->SetTopology((Short_t) Trigger_topology[i]);
1116 HeaderTrig[i]->SetMultiplicity((Short_t) Trigger_multiplicity[i]);
1117 HeaderTrig[i]->SetThreshold(qThreshold[i]);
1118 HeaderTrig[i]->SetAmplitud(Trigger_response_ampl);
1119 HeaderTrig[i]->SetFwhm(Trigger_response_fwhm);
1120 HeaderTrig[i]->SetOverlap(Trigger_overlaping_time);
1121 HeaderTrig[i]->SetGate(Trigger_gate_length);
1122 HeaderTrig[i]->SetElecNoise(Trigger_noise);
1123 }
1124 }
1125 if(Trigger_Loop && Write_McTrig){
1126
1127 int iconcount;
1128 for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++,fthrescount+=Trigger_loop_sthres){
1129 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
1130 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
1131 HeaderTrig[iconcount]->SetTopology((Short_t) isorttopo[itopocount+Trigger_loop_ltop]);
1132 HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
1133 for(int i=0;i<ct_NPixels;i++){
1134 fpixelthres[i]=
1135 ((Float_t)(fthrescount)>=qThreshold[0][i])?
1136 (Float_t)(fthrescount):qThreshold[0][i];
1137 }
1138 HeaderTrig[iconcount]->SetThreshold( fpixelthres);
1139 HeaderTrig[iconcount]->SetAmplitud(Trigger_response_ampl);
1140 HeaderTrig[iconcount]->SetFwhm(Trigger_response_fwhm);
1141 HeaderTrig[iconcount]->SetOverlap(Trigger_overlaping_time);
1142 HeaderTrig[iconcount]->SetGate(Trigger_gate_length);
1143 HeaderTrig[iconcount]->SetElecNoise(Trigger_noise);
1144 iconcount++;
1145 }
1146 }
1147 }
1148 }
1149
1150 // Fill branches for MMcFadcHeader
1151
1152 for(int i=0;i<ct_NPixels;i++){
1153 fadc_elecnoise[i]=FADC_noise;
1154 fadc_diginoise[i]=DIGITAL_noise;
1155 }
1156
1157 Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
1158
1159 if(!Trigger_Loop && Write_McFADC && ct_Number==1){
1160
1161 HeaderFadc[0]->SetShape(0.0);
1162 HeaderFadc[0]->SetAmplitud(FADC_response_ampl*2.35/
1163 sqrt(2*M_PI*FADC_response_fwhm
1164 *FADC_response_fwhm),
1165 FADC_resp_ampl_out*2.35/
1166 sqrt(2*M_PI*FADC_resp_fwhm_out
1167 *FADC_resp_fwhm_out));
1168 HeaderFadc[0]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
1169 HeaderFadc[0]->SetLow2High(FADC_high2low);
1170 HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],
1171 ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
1172 HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
1173 ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
1174 }
1175 if(!Trigger_Loop && Write_McFADC && ct_Number>1){
1176 for(int i=0;i<ct_Number;i++){
1177 Fadc_CT[i]->GetPedestals(&fadc_pedestals[0]);
1178 HeaderFadc[i]->SetShape(0.0);
1179 HeaderFadc[i]->SetAmplitud(FADC_response_ampl*2.35/
1180 sqrt(2*M_PI*FADC_response_fwhm
1181 *FADC_response_fwhm),
1182 FADC_resp_ampl_out*2.35/
1183 sqrt(2*M_PI*FADC_resp_fwhm_out
1184 *FADC_resp_fwhm_out));
1185 HeaderFadc[i]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
1186 HeaderFadc[i]->SetLow2High(FADC_high2low);
1187 HeaderFadc[i]->SetPedestal(&fadc_pedestals[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
1188 HeaderFadc[i]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
1189 ((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
1190 }
1191 }
1192 if(Trigger_Loop && Write_McFADC){
1193 int iconcount;
1194 for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1195 for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
1196 for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
1197 Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
1198 HeaderFadc[iconcount]->SetShape(0.0);
1199 HeaderFadc[iconcount]->SetAmplitud(FADC_response_ampl*2.35/
1200 sqrt(2*M_PI*FADC_response_fwhm
1201 *FADC_response_fwhm),
1202 FADC_resp_ampl_out*2.35/
1203 sqrt(2*M_PI*FADC_resp_fwhm_out
1204 *FADC_resp_fwhm_out));
1205 HeaderFadc[iconcount]->SetFwhm(FADC_response_fwhm,
1206 FADC_resp_fwhm_out);
1207 HeaderFadc[iconcount]->SetLow2High(FADC_high2low);
1208 HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0], ct_NPixels);
1209 HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],
1210 &fadc_diginoise[0],ct_NPixels);
1211 iconcount++;
1212 }
1213 }
1214 }
1215 }
1216
1217 // UnSet flag in pixel 0 (not used for trigger) that indicates if secure pixel
1218 // is active once it has been stored
1219
1220 for (int ict=0; ict<ct_Number;ict++){
1221 qThreshold[ict][0]=999999.99;
1222 }
1223
1224
1225 // create a Tree for the Event data stream
1226 TTree EvtTree("Events","Normal Triggered Events");
1227
1228 if (Write_McEvt && ct_Number==1){
1229
1230 EvtTree.Branch("MMcEvt","MMcEvt",
1231 &McEvt[0]);
1232 }
1233
1234 if (Write_McEvt && ct_Number!=1){
1235 ibr=0;
1236 for(char branchname[10];ibr<numberBranches;ibr++){
1237
1238 sprintf(help,"%i",ibr+1);
1239 strcpy (branchname, "MMcEvt;");
1240 strcat (branchname, & help[0]);
1241 strcat (branchname, ".");
1242 EvtTree.Branch(branchname,"MMcEvt",
1243 &McEvt[ibr]);
1244 }
1245 }
1246
1247 if(!Trigger_Loop && ct_Number==1){
1248
1249 if (Write_RawEvt){
1250 EvtTree.Branch("MRawEvtHeader","MRawEvtHeader",
1251 &EvtHeader[0]);
1252 EvtTree.Branch("MRawEvtData","MRawEvtData",
1253 &EvtData[0]);
1254 }
1255 if (Write_McTrig){
1256 EvtTree.Branch("MMcTrig","MMcTrig",
1257 &McTrig[0]);
1258 }
1259 }
1260 else{
1261
1262 if (Write_McTrig){
1263 ibr=0;
1264 for(char branchname[10];ibr<numberBranches;ibr++){
1265
1266 sprintf(help,"%i",ibr+1);
1267 strcpy (branchname, "MMcTrig;");
1268 strcat (branchname, & help[0]);
1269 strcat (branchname, ".");
1270 EvtTree.Branch(branchname,"MMcTrig",
1271 &McTrig[ibr]);
1272 }
1273 }
1274 }
1275
1276 if ((Trigger_Loop || ct_Number>1) && Write_RawEvt){
1277 ibr=0;
1278 for(char branchname[15];ibr<numberBranches;ibr++){
1279
1280 sprintf(help,"%i",ibr+1);
1281 strcpy (branchname, "MRawEvtHeader;");
1282 strcat (branchname, & help[0]);
1283 strcat (branchname, ".");
1284 EvtTree.Branch(branchname,"MRawEvtHeader",
1285 &EvtHeader[ibr]);
1286 }
1287 ibr=0;
1288 for(char branchname[15];ibr<numberBranches;ibr++){
1289
1290 sprintf(help,"%i",ibr+1);
1291 strcpy (branchname, "MRawEvtData;");
1292 strcat (branchname, & help[0]);
1293 strcat (branchname, ".");
1294 EvtTree.Branch(branchname,"MRawEvtData",
1295 &EvtData[ibr]);
1296 }
1297 }
1298
1299
1300 TApplication theAppTrigger("App", &argc, argv);
1301
1302 if(FADC_Scan){
1303 if (gROOT->IsBatch()) {
1304 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1305 // return 1;
1306 }
1307 }
1308 if(FADC_Scan){
1309 //TApplication theAppFadc("App", &argc, argv);
1310
1311 if (gROOT->IsBatch()) {
1312 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
1313 // return 1;
1314 }
1315 }
1316
1317
1318 // prepare the NSB simulation
1319
1320 // Instance of the Mlons class
1321 MLons lons(Trigger_response_ampl, Trigger_response_fwhm,
1322 FADC_response_ampl,FADC_response_fwhm);
1323
1324 lons.SetSeed(Int_t(get_seeds(1)/get_seeds(0))+1);
1325
1326 lons.SetPath(nsbpathname);
1327
1328 // Instance of the Mlons class
1329 MLons lons_outer(Trigger_response_ampl, Trigger_response_fwhm,
1330 FADC_resp_ampl_out,FADC_resp_fwhm_out);
1331
1332 lons_outer.SetSeed(Int_t(get_seeds(1)/get_seeds(0))+2);
1333
1334 lons_outer.SetPath(nsbpath_outer);
1335
1336 if( simulateNSB ){
1337
1338 //
1339 // Calculate the non-diffuse NSB photoelectron rates
1340 //
1341
1342 // FIXME --- star NSB different for each camera?
1343 log(SIGNATURE,"Produce NSB rates from Star Field");
1344
1345 k = produce_nsbrates( starfieldname,
1346 ((MGeomCam*)(camgeom.UncheckedAt(0))),
1347 nsbrate_phepns,
1348 0);
1349
1350 if (k != 0){
1351 cout << "Error when reading starfield... \nExiting.\n";
1352 exit(1);
1353 }
1354
1355 // calculate diffuse rate correcting for the pixel size and telescope
1356
1357 float factorNSB_ct;
1358 for(int ict=0;ict<ct_Number;ict++){
1359
1360 // First we set the factor due to the mirror size respect to the MAGIC mirror.
1361
1362 switch(GeometryCamera[ict]){
1363 case 1:
1364 case 2:
1365 case 3:
1366 case 4:
1367 factorNSB_ct=1.0;
1368 break;
1369 case 5:
1370 case 6:
1371 case 7:
1372 factorNSB_ct=1000.0/239.0;
1373 break;
1374 case 8:
1375 case 9:
1376 default:
1377 factorNSB_ct=1.0;
1378 break;
1379 }
1380
1381 factorNSB_ct=factorNSB_ct/
1382 ((*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
1383 (*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
1384 PIXEL_SIZE*PIXEL_SIZE);
1385
1386 for(UInt_t ui=0;
1387 ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); ui++){
1388 const Float_t size=
1389 (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD();
1390 diffnsb_phepns[ict][ui] = (Int_t(meanNSB*factorNSB_ct*100*size*size+0.5))/(100.0)*factorqe_NSB[ict];
1391 }
1392 }
1393
1394 // calculate nsb rate including diffuse and starlight
1395 // we also include the extinction effect
1396 for(int ict=0;ict<ct_Number;ict++){
1397 for(j=0;j<iNUMWAVEBANDS;j++){
1398 // calculate the effect of the atmospheric extinction
1399
1400 zenfactor = pow(10., -0.4 * ext[j] );
1401
1402 for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1403 nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
1404 nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
1405 }
1406 }
1407 }
1408 }
1409
1410
1411 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1412 //
1413 // Now 500 empty event with the conditin in which the camera is run
1414 // are simulated. IN this way one gets an estimation of the
1415 // sigma for the pedestal of each FADC channel.
1416 // This sigma computtaion is done assuming any noise taht effects
1417 // the FADC but there is no rotation of the Star Field (otherwise it
1418 // should be done for each event)
1419 //
1420 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1421
1422 for(int ict=0;ict<ct_Number;ict++){
1423 for(UInt_t ui=0;
1424 ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1425 ui++){
1426 fadc_sigma[ui]=0.0;
1427 fadc_sigma_low[ui]=0.0;
1428 }
1429 for(int ie=0;ie<500;ie++){
1430 Fadc_CT[ict]->Reset() ;
1431 if (addElecNoise){
1432 Fadc_CT[ict]->ElecNoise(FADC_noise) ;
1433 }
1434 if(simulateNSB){
1435 for(UInt_t ui=0;
1436 ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1437 ui++){
1438 if(nsb_phepns[ict][ui]>0.0){
1439 if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD()>(*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD()){
1440 k=lons_outer.GetResponse(nsb_phepns[ict][ui],0.01,
1441 & nsb_trigresp[0],
1442 & nsb_fadcresp[0]);
1443 }
1444 else{
1445 k=lons.GetResponse(nsb_phepns[ict][ui],0.01,
1446 & nsb_trigresp[0],& nsb_fadcresp[0]);
1447 }
1448 if(k==0){
1449 cout << "Exiting.\n";
1450 exit(1);
1451 }
1452 Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
1453 }
1454 }
1455 }
1456 Fadc_CT[ict]->Pedestals();
1457 for(UInt_t ui=0;
1458 ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
1459 ui++){
1460 fadc_sigma[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,1))/500.0;
1461 fadc_sigma_low[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,0))/500.0;
1462 }
1463 }
1464 HeaderFadc[ict]->SetPedestalSigma(&fadc_sigma_low[0],&fadc_sigma[0],
1465 ((MGeomCam*)(camgeom.UncheckedAt(ict)))
1466 ->GetNumPixels());
1467 }
1468
1469 //
1470 // Read the reflector file with the Cherenkov data
1471 //
1472
1473 // select input file
1474
1475 if ( Data_From_STDIN ) {
1476
1477 inputfile[0] = stdin;
1478
1479 }
1480 else{
1481
1482 for(int ict=0;ict<ct_Number;ict++){
1483
1484 log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
1485 inputfile[ict] = fopen( inname_CT[ict], "r" );
1486
1487 if ( inputfile[ict] == NULL )
1488 error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
1489 }
1490
1491 }
1492
1493 // get signature, and check it
1494
1495 for(int ict=0;ict<ct_Number;ict++){
1496 if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE){
1497 exit(1);
1498 }
1499 }
1500
1501 // open data file
1502
1503 log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
1504 datafile.open( datname );
1505
1506 if ( datafile.bad() )
1507 error( SIGNATURE, "Cannot open data file: %s\n", datname );
1508
1509 // initializes flag
1510
1511 strcpy( flag, " \0" );
1512 strcpy( flag_new, " \0" );
1513
1514 // allocate space for PMTs numbers of pixels
1515
1516 fnpix = new float [ct_NPixels];
1517
1518
1519 //!@' @#### Main loop.
1520 //@'
1521
1522 // get flag
1523
1524 for(int ict=0;ict<ct_Number;ict++)
1525 fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1526
1527
1528 // loop over the file
1529
1530 still_in_loop = TRUE;
1531
1532 // FIXME --- check if not eof for all input reflector files
1533
1534 while (
1535 ((! Data_From_STDIN) && ( !feof(inputfile[0]) ))
1536 ||
1537 (Data_From_STDIN && still_in_loop)
1538 ) {
1539
1540 // reading .rfl files
1541 if(!isA( flag, FLAG_START_OF_RUN )){
1542
1543 // We break the main loop
1544 cout<<"Warning: Expected start of run flag, but found:"<<flag<<endl;
1545 cout<<" We break the main loop"<<endl;
1546 break;
1547 }
1548 else { // found start of run
1549
1550 for(int ict=0;ict<ct_Number;ict++) {
1551
1552 fread( flag_new, 4, 1, inputfile[ict] );
1553
1554 if(!isA( flag_new, FLAG_START_OF_HEADER)){
1555
1556 // We break the main loop
1557 cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
1558 cout<<" We break the main loop"<<endl;
1559 break;
1560 }
1561
1562 fread((char*)&mcrunh,(SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile[ict] );
1563 }
1564
1565 nshow=0;
1566
1567 for(int ict=0;ict<ct_Number;ict++){
1568 fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
1569 }
1570 while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
1571 for(int ict=0;ict<ct_Number;ict++) {
1572 fread( flag_new, 4, 1, inputfile[ict] );
1573 if(!isA( flag_new, FLAG_EVENT_HEADER)){
1574
1575 // We break while events loop
1576 cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<" "<<ict<<endl;
1577 cout<<" We break the while events loop"<<endl;
1578 break;
1579 }
1580 }
1581
1582
1583 //
1584 // Clear Trigger and Fadc
1585 //
1586 for(int ict=0;ict<ct_Number;ict++) {
1587 Trigger_CT[ict]->Reset() ;
1588 Trigger_CT[ict]->ClearFirst();
1589 Trigger_CT[ict]->ClearZero();
1590 Fadc_CT[ict]->Reset() ;
1591 }
1592
1593 ++nshow;
1594 if((nshow+ntshow+1)%100 == 1)
1595 log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
1596
1597 // get MCEventHeader
1598 if (reflector_file_version<6){
1599 for(int ict=0;ict<ct_Number;ict++)
1600 fread( (char*)&mcevth[ict], mcevth[ict].mysize(),
1601 1, inputfile[ict] );
1602 }
1603 else{
1604 for(int ict=0;ict<ct_Number;ict++)
1605 fread( (char*)&mcevth_2[ict], mcevth_2[ict].mysize(),
1606 1, inputfile[ict] );
1607 }
1608
1609 // calculate impact parameter (shortest distance betwee the original
1610 // trajectory of the primary (assumed shower-axis) and the
1611 // direction where the telescope points to
1612 //
1613 // we use the following equation, given that the shower core position
1614 // is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
1615 // and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
1616 // and (l2,m2,n2)
1617 //
1618 // | |
1619 // | x1-x2 y1-y2 z1-z2 |
1620 // | |
1621 // + | l1 m1 n1 |
1622 // - | |
1623 // | l2 m2 n2 |
1624 // | |
1625 // dist = ------------------------------------ ( > 0 )
1626 // [ |l1 m1|2 |m1 n1|2 |n1 l1|2 ]1/2
1627 // [ | | + | | + | | ]
1628 // [ |l2 m2| |m2 n2| |n2 l2| ]
1629 //
1630 // playing a little bit, we get this reduced for in our case:
1631 //
1632 //
1633 // dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
1634 // [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
1635 // 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
1636
1637 // read the direction of the incoming shower
1638 // It is done only for one telescope since it is suposed
1639 // to be the same shower for all of them
1640 if (reflector_file_version<6){
1641 thetashw = mcevth[0].get_theta();
1642 phishw = mcevth[0].get_phi();
1643 }
1644 else{
1645 thetashw = mcevth_2[0].get_theta();
1646 phishw = mcevth_2[0].get_phi();
1647 }
1648
1649 // calculate vector for shower
1650
1651 l1 = sin(thetashw)*cos(phishw);
1652 m1 = sin(thetashw)*sin(phishw);
1653 n1 = cos(thetashw);
1654
1655 if (reflector_file_version<6){
1656
1657 mcevth[0].get_core(&coreX, &coreY);
1658
1659 for(int ict=0;ict<ct_Number;ict++){
1660 // read the deviation of the telescope with respect to the shower
1661 mcevth[ict].get_deviations ( &thetaCT[ict], &phiCT[ict] );
1662
1663 if ( (thetaCT[ict] == 0.) && (phiCT[ict] == 0.) ) {
1664
1665 // CT was looking to the source (both lines are parallel)
1666 // therefore, we calculate the impact parameter as the distance
1667 // between the CT axis and the core position
1668
1669 impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1670 thetaCT[ict] += thetashw;
1671 phiCT[ict] += phishw;
1672 } else {
1673
1674 // the shower comes off-axis
1675
1676 // obtain with this the final direction of the CT
1677
1678 thetaCT[ict] += thetashw;
1679 phiCT[ict] += phishw;
1680
1681 // calculate vector for telescope
1682
1683 l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
1684 m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
1685 n2 = cos(thetaCT[ict]);
1686
1687 num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
1688 den = (SQR(l1*m2 - l2*m1) +
1689 SQR(m1*n2 - m2*n1) +
1690 SQR(n1*l2 - n2*l1));
1691 den = sqrt(den);
1692 impactD[ict] = fabs(num)/den;
1693 }
1694 }
1695 }
1696 else{
1697 mcevth_2[0].get_core(&coreX, &coreY);
1698 for(int ict=0;ict<ct_Number;ict++){
1699 thetaCT[ict]=mcevth_2[ict].get_theta_CT();
1700 phiCT[ict]=mcevth_2[ict].get_phi_CT();
1701 if ( (thetaCT[ict] == thetashw) && (phiCT[ict] == phishw) ) {
1702
1703 // CT was looking to the source (both lines are parallel)
1704 // therefore, we calculate the impact parameter as the distance
1705 // between the CT axis and the core position
1706
1707 impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
1708
1709 } else {
1710
1711 // the shower comes off-axis
1712
1713 // calculate vector for telescope
1714
1715 l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
1716 m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
1717 n2 = cos(thetaCT[ict]);
1718
1719 num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
1720 den = (SQR(l1*m2 - l2*m1) +
1721 SQR(m1*n2 - m2*n1) +
1722 SQR(n1*l2 - n2*l1));
1723 den = sqrt(den);
1724
1725 impactD[ict] = fabs(num)/den;
1726
1727 }
1728 }
1729 }
1730
1731 // energy cut
1732
1733 if ( Select_Energy ) {
1734 if (reflector_file_version<6)
1735 if (( mcevth[0].get_energy() < Select_Energy_le ) ||
1736 ( mcevth[0].get_energy() > Select_Energy_ue )) {
1737 log(SIGNATURE, "select_energy: shower rejected.\n");
1738 continue;
1739 }
1740 else
1741 if (( mcevth_2[0].get_energy() < Select_Energy_le ) ||
1742 ( mcevth_2[0].get_energy() > Select_Energy_ue )) {
1743 log(SIGNATURE, "select_energy: shower rejected.\n");
1744 continue;
1745 }
1746 }
1747
1748 inumphe=0;
1749
1750 for(int ict=0;ict<ct_Number;ict++){
1751
1752 // Read first and last time and put inumphe_CT[0] to 0
1753
1754 if (reflector_file_version<6)
1755 mcevth[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1756 else
1757 mcevth_2[ict].get_times(&arrtmin_ns,&arrtmax_ns);
1758
1759 inumphe_CT[ict]=0;
1760
1761 // read the photons and produce the photoelectrons
1762
1763 k = produce_phes( inputfile[ict],
1764 ((MGeomCam*)(camgeom.UncheckedAt(ict))),
1765 WAVEBANDBOUND1,
1766 WAVEBANDBOUND6,
1767 Trigger_CT[ict], // will be changed by the function!
1768 Fadc_CT[ict], // will be changed by the function!
1769 &inumphe_CT[ict], // important for later: the size of photoe[]
1770 fnpix, // will be changed by the function!
1771 &ncph[ict], // will be changed by the function!
1772 &arrtmin_ns, // will be changed by the function!
1773 &arrtmax_ns, // will be changed by the function!
1774 ict);
1775
1776 inumphe=(inumphe<inumphe_CT[ict])?inumphe_CT[ict]:inumphe;
1777
1778 if( k != 0 ){ // non-zero returnvalue means error
1779 cout << "Exiting.\n";
1780 exit(1);
1781 }
1782 }
1783
1784 // NSB simulation
1785
1786 if(simulateNSB && nphe2NSB<=inumphe){
1787
1788 if(Starfield_rotate){
1789
1790 // Introduction rho angle
1791
1792 zenith = thetashw;
1793 azimutal = phishw;
1794 C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
1795 C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
1796 C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
1797 rho = acos( C1/C2 );
1798
1799 if ( sin(azimutal) < 0)
1800 {
1801 rho = 2 * 3.14159 - rho;
1802 }
1803 else
1804 {
1805 rho = rho;
1806 }
1807 rho = rho*180/3.14159;
1808
1809 // Rotation of the NSB
1810 // FIXME --- We should rotate for all cameras. Is it always the same rho?
1811 for(int ict=0;ict<ct_Number;ict++){
1812 k = size_rotated(
1813 &nsb_phepns_rotated[ict][0],
1814 nsb_phepns[ict],
1815 rho);
1816 }
1817 }
1818
1819 // Fill trigger and fadc response in the trigger class from the database
1820 for(int ict=0;ict<ct_Number;ict++){
1821
1822 for(UInt_t ui=0;ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1823 if(nsb_phepns_rotated[ict][ui]>0.0){
1824 if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD()>(*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD()){
1825 k=lons_outer.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
1826 & nsb_trigresp[0],
1827 & nsb_fadcresp[0]);
1828 }
1829 else{
1830 k=lons.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
1831 & nsb_trigresp[0],& nsb_fadcresp[0]);
1832 }
1833 if(k==0){
1834 cout << "Exiting.\n";
1835 exit(1);
1836 }
1837 Trigger_CT[ict]->AddNSB(ui,nsb_trigresp);
1838 Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
1839 }
1840 }
1841 }
1842
1843 }// end if(simulateNSB && nphe2NSB<=inumphe_CT[0]) ...
1844
1845 for(int ict=0;ict<ct_Number;ict++){
1846 inumphensb[ict]=0;
1847 for (UInt_t ui=0;ui<
1848 ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
1849 inumphensb[ict]+=nsb_phepns[ict][ui]*TOTAL_TRIGGER_TIME;
1850 }
1851 ntcph[ict]+=ncph[ict];
1852 if ((nshow+ntshow+1)%100 == 1){
1853 log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1854 ncph[ict], ntcph[ict]);
1855
1856 cout << "Total number of phes in CT "<<ict<<": "
1857 << inumphe_CT[ict]<<" (+ ";
1858 cout<<inumphensb[ict]<<" mean expected number from NSB)"<<endl;
1859 }
1860 }
1861
1862 // skip it ?
1863
1864 int i;
1865 for (i=0; i<nSkip; ++i ) {
1866 if (Skip[i] == (nshow+ntshow)) {
1867 i = -1;
1868 break;
1869 }
1870 }
1871
1872 // if after the previous loop, the exit value of i is -1
1873 // then the shower number is in the list of showers to be
1874 // skipped
1875
1876 if (i == -1) {
1877 log(SIGNATURE, "\t\tskipped!\n");
1878 continue;
1879 }
1880
1881 //++++++++++++++++++++++++++++++++++++++++++++++++++
1882 // at this point we have a camera full of
1883 // ph.e.s
1884 // we should first apply the trigger condition,
1885 // and if there's trigger, then clean the image,
1886 // calculate the islands statistics and the
1887 // other parameters of the image (Hillas' parameters
1888 // and so on).
1889 //--------------------------------------------------
1890
1891 // TRIGGER HERE
1892
1893 // We should simulate the AC coupling behaviour:
1894 // For the FADC it is only done for the NSB while producinr
1895 // the StarResponse database.
1896 // For the trigger is done in the Trigger.Diskriminate(), which
1897 // is called later (it should be seperated to speed up the program
1898 //
1899
1900 // now the noise of the electronic
1901 // (preamps, optical transmission,..) is introduced.
1902 // This is done inside the class MTrigger by the method ElecNoise.
1903 //
1904
1905 for(int ict=0;ict<ct_Number;ict++) {
1906 if (addElecNoise && nphe2NSB<=inumphe){
1907
1908 Trigger_CT[ict]->ElecNoise(Trigger_noise) ;
1909
1910 Fadc_CT[ict]->ElecNoise(FADC_noise) ;
1911 }
1912 }
1913
1914 // now a shift in the fadc signal due to the pedestlas is
1915 // intorduced
1916 // This is done inside the class MFadc by the method Pedestals
1917 for(int ict=0;ict<ct_Number;ict++) {
1918 Fadc_CT[ict]->Pedestals();
1919 }
1920
1921 // We study several trigger conditons
1922 if(Trigger_Loop){
1923
1924 // Set to zero the flag to know if some conditon has triggered
1925 btrigger=0;
1926 flagstoring = 0;
1927
1928 // Loop over trigger threshold
1929 int iconcount;
1930 for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
1931 for (int i=0;i<ct_NPixels;i++){
1932 fpixelthres[i]=
1933 ((Float_t)(fthrescount)>=qThreshold[0][i])?
1934 (Float_t)(fthrescount):qThreshold[0][i];
1935
1936 // Rise the discrimnator threshold to avoid huge rates
1937
1938 if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
1939 for(int ii=0;ii<ct_NPixels;ii++)
1940 if( nsb_phepns_rotated[0][ii]>riseDiskThres)
1941 fpixelthres[ii]=secureDiskThres;
1942 }
1943 Trigger_CT[0]->SetThreshold(fpixelthres);
1944
1945 Trigger_CT[0]->Diskriminate();
1946
1947
1948 //
1949 // look if in all the signals in the trigger signal branch
1950 // is a possible Trigger. Therefore we habe to diskriminate all
1951 // the simulated analog signals (Method Diskriminate in class
1952 // MTrigger). We look simultanously for the moments at which
1953 // there are more than TRIGGER_MULTI pixels above the
1954 // CHANNEL_THRESHOLD.
1955 //
1956
1957 // Set trigger flags to zero
1958 Lev0=0;
1959 Lev1=0;
1960
1961 // loop over multiplicity of trigger configuration
1962 for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
1963 Trigger_CT[0]->SetMultiplicity(imulticount);
1964 Trigger_CT[0]->ClearZero();
1965
1966 Lev0=(Short_t) Trigger_CT[0]->ZeroLevel();
1967 if (Lev0>0 || Write_All_Images || btrigger){
1968
1969 // loop over topologies
1970 for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
1971 Lev1=0;
1972
1973 if(itopocount==0 && imulticount>7) continue;
1974 //COBB if(itopocount==2 && imulticount<3) continue;
1975 // It only makes to look for a different topology
1976 // if there are 3 or more N pixels.
1977 if(imulticount<3)
1978 Trigger_CT[0]->SetTopology(1);
1979 else
1980 {
1981 // We should be careful that topologies are sort from
1982 // the less to the more restrictive one.
1983 Trigger_CT[0]->SetTopology(isorttopo[itopocount]);
1984 }
1985 Trigger_CT[0]->ClearFirst();
1986
1987 //
1988 // Start the First Level Trigger simulation
1989 //
1990 if(Lev0!=0)
1991 Lev1=Trigger_CT[0]->FirstLevel();
1992 if(Lev1>0) {
1993 btrigger= 1;
1994 ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++;
1995 }
1996
1997 Lev0=1;
1998 Int_t NumImages = Lev1;
1999 if(Lev1==0 && (Write_All_Images || btrigger)){
2000 btrigger= 1;
2001 NumImages=1;
2002 Lev0=0;
2003 }
2004
2005 for (Int_t ii=0;ii<NumImages;ii++){
2006 if (Write_McTrig){
2007 McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
2008 McTrig[iconcount]->SetTime(Trigger_CT[0]->GetFirstLevelTime(ii),ii+1);
2009 Trigger_CT[0]->GetMapDiskriminator(trigger_map);
2010 McTrig[iconcount]->SetMapPixels(trigger_map,ii);
2011 }
2012 //
2013 // fill inside the class fadc the member output
2014 //
2015
2016 Fadc_CT[0]->TriggeredFadc(Trigger_CT[0]->GetFirstLevelTime(ii));
2017
2018 if( Write_RawEvt ){
2019 //
2020 // Fill the header of this event
2021 //
2022
2023 EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
2024 // fill pixel information
2025 if (Lev1 || Write_All_Images){
2026 if (addElecNoise) Fadc_CT[0]->DigitalNoise();
2027 for(UInt_t i=0;
2028 i<((MGeomCam*)(camgeom.UncheckedAt(0)))
2029 ->GetNumPixels();i++){
2030//
2031// AM 15 01 2004: commented out "continue" statement, so that also pixels
2032// with no C-photons will be written to the output in case the camera
2033// is run with no noise.
2034// if(!Fadc_CT[0]->IsPixelUsed(i)) continue;
2035//
2036 for (j=0;j<FADC_SLICES;j++){
2037 fadcValues->AddAt(Fadc_CT[0]->GetFadcSignal(i,j),j);
2038 fadcValuesLow->AddAt(Fadc_CT[0]->GetFadcLowGainSignal(i,j),j);
2039 }
2040 EvtData[iconcount]->AddPixel(i,fadcValues,0);
2041 EvtData[iconcount]->AddPixel(i,fadcValuesLow,kTRUE);
2042 }
2043 }
2044 }
2045 }
2046 //
2047 // Increase counter of analised trigger conditions
2048 //
2049 iconcount++;
2050 }
2051 }
2052 else{
2053 break;
2054 }
2055 }
2056 if (!btrigger) break;
2057 }
2058 if (btrigger){
2059
2060 //
2061 // fill the MMcEvt with all information
2062 //
2063
2064 if(!flagstoring)
2065 nstoredevents++;
2066 flagstoring = 1;
2067
2068 if (Write_McEvt) {
2069 Float_t ftime, ltime;
2070 if (reflector_file_version<6){
2071 mcevth[0].get_times(&ftime, &ltime);
2072 McEvt[0]->Fill( 0,
2073 (UShort_t) mcevth[0].get_primary() ,
2074 mcevth[0].get_energy(),
2075 -1.0,
2076 -1.0,
2077 -1.0,
2078 mcevth[0].get_theta(),
2079 mcevth[0].get_phi(),
2080 mcevth[0].get_core(),
2081 mcevth[0].get_coreX(),
2082 mcevth[0].get_coreY(),
2083 impactD[0],
2084 phiCT[0],
2085 thetaCT[0],
2086 ftime,
2087 ltime,
2088 0,
2089 0,
2090 0,
2091 0,
2092 0,
2093 0,
2094 0,
2095 (UInt_t)mcevth[0].get_CORSIKA(),
2096 (UInt_t)mcevth[0].get_AtmAbs(),
2097 (UInt_t)(mcevth[0].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()),
2098 (UInt_t) ncph[0],
2099 (UInt_t) inumphe_CT[0],
2100 (UInt_t) inumphensb[0]+inumphe_CT[0],
2101 -1.0,
2102 -1.0,
2103 -1.0);
2104 }
2105 else{
2106 Float_t Nmax, t0, tmax, a, b, c, chi2;
2107 mcevth_2[0].get_times(&ftime, &ltime);
2108 chi2=mcevth_2[0].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2109 McEvt[0]->Fill((UInt_t) mcevth_2[0].get_evt_number(),
2110 (UShort_t) mcevth_2[0].get_primary() ,
2111 mcevth_2[0].get_energy(),
2112 mcevth_2[0].get_thick0(),
2113 mcevth_2[0].get_first_target(),
2114 mcevth_2[0].get_z_first_int(),
2115 mcevth_2[0].get_theta(),
2116 mcevth_2[0].get_phi(),
2117 mcevth_2[0].get_core(),
2118 mcevth_2[0].get_coreX(),
2119 mcevth_2[0].get_coreY(),
2120 impactD[0],
2121 mcevth_2[0].get_phi_CT(),
2122 mcevth_2[0].get_theta_CT(),
2123 ftime,
2124 ltime,
2125 Nmax,
2126 t0,
2127 tmax,
2128 a,
2129 b,
2130 c,
2131 chi2,
2132 (UInt_t)mcevth_2[0].get_CORSIKA(),
2133 (UInt_t)mcevth_2[0].get_AtmAbs(),
2134 (UInt_t)(mcevth_2[0].get_MirrAbs()+mcevth_2[0].get_OutOfMirr()+mcevth_2[0].get_BlackSpot()),
2135 (UInt_t) ncph[0],
2136 (UInt_t) inumphe_CT[0],
2137 (UInt_t) inumphensb[0]+inumphe_CT[0],
2138 mcevth_2[0].get_ElecFraction(),
2139 mcevth_2[0].get_MuonFraction(),
2140 mcevth_2[0].get_OtherFraction());
2141 }
2142 }
2143 // Fill the Tree with the current leaves of each branch
2144 i=EvtTree.Fill() ;
2145
2146 // Clear the branches
2147 if(Write_McTrig){
2148 for(int i=0;i<numberBranches;i++){
2149 McTrig[i]->Clear() ;
2150 }
2151 }
2152 if( Write_RawEvt ){
2153 for(int i=0;i<numberBranches;i++){
2154 EvtHeader[i]->Clear() ;
2155 EvtData[i]->DeletePixels();
2156 }
2157 }
2158 if (Write_McEvt)
2159 McEvt[0]->Clear() ;
2160 }
2161 }
2162 // We study a single trigger condition
2163 else {
2164
2165 // Set to zero the flag to know if some conditon has triggered
2166 btrigger=0;
2167 flagstoring = 0;
2168
2169 for(int ict=0;ict<ct_Number;ict++){
2170
2171 // Setting trigger conditions
2172 Trigger_CT[ict]->SetMultiplicity(Trigger_multiplicity[ict]);
2173 Trigger_CT[ict]->SetTopology(Trigger_topology[ict]);
2174 for (int i=0;i<ct_NPixels;i++)
2175 fpixelthres[i]=qThreshold[ict][i];
2176
2177 // Rise the discrimnator threshold to avoid huge rates
2178 if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
2179 for(int ii=0;ii<ct_NPixels;ii++){
2180 if( nsb_phepns_rotated[ict][ii]>riseDiskThres)
2181 fpixelthres[ii]=secureDiskThres;
2182 }
2183
2184 Trigger_CT[ict]->SetThreshold(fpixelthres);
2185
2186 Trigger_CT[ict]->Diskriminate() ;
2187
2188 //
2189 // look if in all the signals in the trigger signal branch
2190 // is a possible Trigger. Therefore we habe to diskriminate all
2191 // the simulated analog signals (Method Diskriminate in class
2192 // MTrigger). We look simultanously for the moments at which
2193 // there are more than TRIGGER_MULTI pixels above the
2194 // CHANNEL_THRESHOLD.
2195 //
2196
2197 Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ;
2198
2199 Lev1MT[ict] = 0 ;
2200
2201 //
2202 // Start the First Level Trigger simulation
2203 //
2204
2205 if ( Lev0MT[ict] > 0 || Write_All_Images) {
2206
2207 Lev1MT[ict]= Trigger_CT[ict]->FirstLevel();
2208 }
2209 if (Lev1MT[ict]>0){
2210 ++ntrigger[ict];
2211 }
2212 }
2213
2214 Int_t NumImages = 0;
2215 Int_t CT_triggered=0;
2216 for(int ict=0;ict<ct_Number;ict++){
2217 if(NumImages==0 && Lev1MT[ict]>0)
2218 CT_triggered=ict;
2219 NumImages = (NumImages>=Lev1MT[ict])?NumImages:1;
2220 Lev0MT[ict]=1;
2221 if (Lev1MT[ict]==0 && Write_All_Images){
2222 NumImages=1;
2223 Lev0MT[ict]=0;
2224 }
2225 }
2226
2227 for(int ict=0;ict<ct_Number;ict++){
2228 for(Int_t ii=0;ii<NumImages;ii++){
2229
2230 btrigger=1;
2231
2232 // Loop over different level one triggers
2233
2234 //
2235 // fill inside class fadc the member output
2236 //
2237 if(Lev1MT[ict]>0)
2238 Fadc_CT[ict]->
2239 TriggeredFadc(Trigger_CT[ict]->GetFirstLevelTime(ii));
2240 else
2241 Fadc_CT[ict]->
2242 TriggeredFadc(Trigger_CT[CT_triggered]->GetFirstLevelTime(ii));
2243
2244 if(!flagstoring)
2245 nstoredevents++;
2246 flagstoring = 1;
2247
2248 if (Write_McTrig){
2249 McTrig[ict]->SetFirstLevel (Lev1MT[ict]);
2250 McTrig[ict]
2251 ->SetTime(Trigger_CT[ict]->GetFirstLevelTime(ii),ii+1);
2252 Trigger_CT[ict]->GetMapDiskriminator(trigger_map);
2253 McTrig[ict]->SetMapPixels(trigger_map,ii);
2254 }
2255
2256 // Fill Evt information
2257
2258 if (Write_RawEvt){
2259
2260 //
2261 // Fill the header of this event
2262 //
2263
2264 EvtHeader[ict]
2265 ->FillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ;
2266
2267 // fill pixel information
2268
2269 if (Lev1MT[ict] || Write_All_Images){
2270 if (addElecNoise) Fadc_CT[ict]->DigitalNoise();
2271 for(UInt_t i=0;
2272 i<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
2273 i++){
2274//
2275// AM 15 01 2004 Commented out "continue" statement, so that also pixels
2276// containing no C-photons will be written to the output in case of running
2277// camera with no noise added to the signal.
2278// if(!Fadc_CT[ict]->IsPixelUsed(i)) continue;
2279//
2280 for (j=0;j<FADC_SLICES;j++){
2281 fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(i,j),j);
2282 fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(i,j),j);
2283 }
2284 EvtData[ict]->AddPixel(i,fadcValues,0);
2285 EvtData[ict]->AddPixel(i,fadcValuesLow,kTRUE);
2286 }
2287 }
2288 }
2289 }
2290 //
2291 // if a first level trigger occurred, then
2292 // 1. do some other stuff (not implemented)
2293 // 2. start the gui tool
2294
2295 if(FADC_Scan){
2296 if ( Lev0MT[ict] > 0 ) {
2297 Fadc_CT[ict]->ShowSignal( McEvt[ict], (Float_t) 60. ) ;
2298 }
2299 }
2300
2301 if(Trigger_Scan){
2302 if ( Lev0MT[ict] > 0 ) {
2303 Trigger_CT[ict]->ShowSignal(McEvt[ict]) ;
2304 }
2305 }
2306
2307 } // end CT loop
2308
2309 // If there is trigger in some telecope or we store all showers
2310 if(btrigger){
2311 if (Write_McEvt){
2312 //
2313 // fill the MMcEvt with all information
2314 //
2315
2316 for (int ict=0;ict<ct_Number;ict++){
2317 Float_t ftime, ltime;
2318 if (reflector_file_version<6){
2319 mcevth[ict].get_times(&ftime, &ltime);
2320 McEvt[ict]->Fill( 0,
2321 (UShort_t) mcevth[0].get_primary() ,
2322 mcevth[ict].get_energy(),
2323 -1.0,
2324 -1.0,
2325 -1.0,
2326 mcevth[ict].get_theta(),
2327 mcevth[ict].get_phi(),
2328 mcevth[ict].get_core(),
2329 mcevth[ict].get_coreX(),
2330 mcevth[ict].get_coreY(),
2331 impactD[ict],
2332 phiCT[ict],
2333 thetaCT[ict],
2334 ftime,
2335 ltime,
2336 0,
2337 0,
2338 0,
2339 0,
2340 0,
2341 0,
2342 0,
2343 (UInt_t)mcevth[ict].get_CORSIKA(),
2344 (UInt_t)mcevth[ict].get_AtmAbs(),
2345 (UInt_t)(mcevth[ict].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()),
2346 (UInt_t) ncph[ict],
2347 (UInt_t) inumphe_CT[ict],
2348 (UInt_t) inumphensb[ict]+inumphe_CT[ict],
2349 -1.0,
2350 -1.0,
2351 -1.0);
2352 }
2353 else{
2354 Float_t Nmax, t0, tmax, a, b, c, chi2;
2355 mcevth_2[ict].get_times(&ftime, &ltime);
2356 chi2=mcevth_2[ict].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
2357 McEvt[ict]->Fill( (UInt_t) mcevth_2[ict].get_evt_number(),
2358 (UShort_t) mcevth_2[ict].get_primary() ,
2359 mcevth_2[ict].get_energy(),
2360 mcevth_2[ict].get_thick0(),
2361 mcevth_2[ict].get_first_target(),
2362 mcevth_2[ict].get_z_first_int(),
2363 mcevth_2[ict].get_theta(),
2364 mcevth_2[ict].get_phi(),
2365 mcevth_2[ict].get_core(),
2366 mcevth_2[ict].get_coreX(),
2367 mcevth_2[ict].get_coreY(),
2368 impactD[ict],
2369 mcevth_2[ict].get_phi_CT(),
2370 mcevth_2[ict].get_theta_CT(),
2371 ftime,
2372 ltime,
2373 Nmax,
2374 t0,
2375 tmax,
2376 a,
2377 b,
2378 c,
2379 chi2,
2380 (UInt_t)mcevth_2[ict].get_CORSIKA(),
2381 (UInt_t)mcevth_2[ict].get_AtmAbs(),
2382 (UInt_t) (mcevth_2[ict].get_MirrAbs()+mcevth_2[ict].get_OutOfMirr()+mcevth_2[ict].get_BlackSpot()),
2383 (UInt_t) ncph[ict],
2384 (UInt_t) inumphe_CT[ict],
2385 (UInt_t) inumphensb[ict]+inumphe_CT[ict],
2386 mcevth_2[ict].get_ElecFraction(),
2387 mcevth_2[ict].get_MuonFraction(),
2388 mcevth_2[ict].get_OtherFraction());
2389 }
2390 }
2391 }
2392 // We don not count photons out of the camera.
2393
2394
2395 //
2396 // write it out to the file outfile
2397 //
2398
2399 EvtTree.Fill() ;
2400 }
2401
2402 // clear all
2403 for(int ict=0;ict<ct_Number;ict++){
2404 if (Write_RawEvt) EvtHeader[ict]->Clear() ;
2405 if (Write_RawEvt) EvtData[ict]->DeletePixels();
2406 if (Write_McTrig) McTrig[ict]->Clear() ;
2407 if (Write_McEvt) McEvt[ict]->Clear() ;
2408 }
2409 }
2410
2411#ifdef __DEBUG__
2412 printf("\n");
2413
2414 for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
2415
2416 for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
2417
2418 if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
2419
2420 if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
2421
2422 printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
2423 (int)pixels[ici][icj][PIXNUM],
2424 pixels[ici][icj][PIXX],
2425 pixels[ici][icj][PIXY],
2426 fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
2427
2428 }
2429 }
2430 }
2431 }
2432
2433 for (int i=0;
2434 i<((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels(); ++i) {
2435 printf("%d (%d): ", i, npixneig[i]);
2436 for (j=0; j<npixneig[i]; ++i)
2437 printf(" %d", pixneig[i][j]);
2438 printf("\n");
2439 }
2440
2441#endif // __DEBUG__
2442
2443
2444 // We search the maximum impact parameter fo the simualted showers
2445 maxpimpact=maxpimpact<impactD[0]?impactD[0]:maxpimpact;
2446
2447 // look for the next event
2448
2449 for(int ict=0;ict<ct_Number;ict++)
2450 fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2451
2452 } // end while there is a next event
2453
2454 if( !isA( flag, FLAG_END_OF_RUN )){
2455 error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
2456 break;
2457 }
2458 else { // found end of run
2459 ntshow += nshow;
2460 log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
2461
2462 //fread( flag, SIZE_OF_FLAGS, 1, inputfile );
2463 for(int ict=0;ict<ct_Number;ict++)
2464 fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2465
2466 if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
2467 log(SIGNATURE, "End of file . . .\n");
2468 still_in_loop = FALSE;
2469 log(SIGNATURE, "Reading ASCII files at the end of the reflector file. . .\n");
2470 for(int ict=0;ict<ct_Number;ict++){
2471 read_ascii(inputfile[ict], McConfigRunHeader[ict]);
2472
2473 if (get_sigma_xy_cm_spot() > 0.)
2474 {
2475 Float_t spotsigma = McConfigRunHeader[ict]->GetPointSpread();
2476 Float_t newsigma = sqrt(spotsigma * spotsigma +
2477 get_sigma_xy_cm_spot() * get_sigma_xy_cm_spot());
2478 McConfigRunHeader[ict]->SetPointSpread(newsigma);
2479 }
2480 }
2481 if ((! Data_From_STDIN) && ( !feof(inputfile[0]) )){
2482
2483 // we have concatenated input files.
2484 // get signature of the next part and check it.
2485
2486 if((reflector_file_version=check_reflector_file( inputfile[0] ))==FALSE){
2487 log(SIGNATURE, "Next file is not recognised as a reflector file.\n");
2488 log(SIGNATURE, "Stopping ...\n");
2489 break;
2490 }
2491
2492 }
2493
2494 for(int ict=0;ict<ct_Number;ict++)
2495 fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
2496
2497 } // end if found end of file
2498 } // end if found end of run
2499
2500 } // end if else found start of run
2501 } // end big while loop
2502
2503 //<@ Finally we should fill th McRunHeader
2504
2505 Float_t heights[10];
2506 time_t ltime;
2507 Float_t ftime;
2508 Float_t rnum;
2509 Float_t viewcone[2]={0,0};
2510
2511 get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
2512 if (reflector_file_version<6){
2513 telestheta=-10.0;
2514 telesphi=-10.0;
2515 mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2516 mcevth[0].get_phi_range(&shphimin,&shphimax);
2517 mcevth[0].get_theta_range(&shthetamin, &shthetamax);
2518 mcevth[0].get_phi_range(&shphimin,&shphimax);
2519 corsika=UInt_t(mcevth[0].get_VersionPGM()*1000);
2520 for (int i=0; i< 10;i++)
2521 heights[i]=mcevth[0].get_HeightLev (i);
2522 rnum=mcevth[0].get_RunNumber();
2523 }
2524 else{
2525 telestheta=mcevth_2[0].get_theta_CT();
2526 telesphi=mcevth_2[0].get_phi_CT();
2527 mcevth_2[0].get_theta_range(&shthetamin, &shthetamax);
2528 mcevth_2[0].get_phi_range(&shphimin,&shphimax);
2529 corsika=UInt_t(mcevth_2[0].get_VersionPGM()*1000);
2530 for (int i=0; i< 10;i++)
2531 heights[i]=mcevth_2[0].get_HeightLev (i);
2532 rnum=mcevth_2[0].get_RunNumber();
2533 mcevth_2[0].get_viewcone(&viewcone[0],&viewcone[1]);
2534 }
2535
2536 get_source_off(&sofftheta,&soffphi);
2537 if(!Trigger_Loop) icontrigger=0;
2538 time (&ltime);
2539 ftime = ((Float_t)ltime)/1000;
2540
2541 if (reflector_file_version<6)
2542 McRunHeader->Fill(rnum,
2543 (UInt_t) 0,
2544 mcevth[0].get_DateRun(),
2545 ftime,
2546 icontrigger,
2547 !Write_All_Images,
2548 Write_McEvt,
2549 Write_McTrig,
2550 Write_McFADC,
2551 Write_RawEvt,
2552 addElecNoise,
2553 ct_NPixels,
2554 (UInt_t)ntshow,
2555 (UInt_t)nstoredevents,
2556 0,
2557 sfRaH,
2558 sfRaM,
2559 sfRaS,
2560 sfDeD,
2561 sfDeM,
2562 sfDeS,
2563 meanNSB,
2564 telestheta,
2565 telesphi,
2566 sofftheta,
2567 soffphi,
2568 shthetamax,
2569 shthetamin,
2570 shphimax,
2571 shphimin,
2572 maxpimpact,
2573 mcevth[0].get_CWaveLower(),
2574 mcevth[0].get_CWaveUpper(),
2575 mcevth[0].get_slope(),
2576 1,
2577 heights,
2578 corsika,
2579 (UInt_t)(reflector_file_version*100),
2580 (UInt_t)(VERSION*100),
2581 0);
2582 else
2583 McRunHeader->Fill(rnum,
2584 (UInt_t) 0,
2585 mcevth_2[0].get_DateRun(),
2586 ftime,
2587 icontrigger,
2588 !Write_All_Images,
2589 Write_McEvt,
2590 Write_McTrig,
2591 Write_McFADC,
2592 Write_RawEvt,
2593 addElecNoise,
2594 ct_NPixels,
2595 (UInt_t)ntshow,
2596 (UInt_t)nstoredevents,
2597 0,
2598 sfRaH,
2599 sfRaM,
2600 sfRaS,
2601 sfDeD,
2602 sfDeM,
2603 sfDeS,
2604 meanNSB,
2605 telestheta,
2606 telesphi,
2607 sofftheta,
2608 soffphi,
2609 shthetamax,
2610 shthetamin,
2611 shphimax,
2612 shphimin,
2613 maxpimpact,
2614 mcevth_2[0].get_CWaveLower(),
2615 mcevth_2[0].get_CWaveUpper(),
2616 mcevth_2[0].get_slope(),
2617 1,
2618 heights,
2619 corsika,
2620 (UInt_t)(reflector_file_version*100),
2621 (UInt_t)(VERSION*100),
2622 0);
2623 // Fill some missing values for MRawRunHeader
2624
2625 RunHeader->SetRunNumber(rnum);
2626 RunHeader->SetNumEvents(nstoredevents);
2627
2628 // Fill MMcCorsikaRunHeader
2629
2630 Float_t constantC[50];
2631 mcrunh.get_constantC(&constantC[0]);
2632 Float_t constantCKA[40];
2633 mcrunh.get_constantCKA(&constantCKA[0]);
2634 Float_t constantCETA[5];
2635 mcrunh.get_constantCETA(&constantCETA[0]);
2636 Float_t constantCSTRBA[11];
2637 mcrunh.get_constantCSTRBA(&constantCSTRBA[0]);
2638 Float_t constantAATM[5];
2639 mcrunh.get_constantAATM(&constantAATM[0]);
2640 Float_t constantBATM[5];
2641 mcrunh.get_constantBATM(&constantBATM[0]);
2642 Float_t constantCATM[5];
2643 mcrunh.get_constantCATM(&constantCATM[0]);
2644 Float_t constantNFL[4];
2645 mcrunh.get_constantNFL(&constantNFL[0]);
2646
2647 if(reflector_file_version>5)
2648 McCorsikaRunHeader->Fill(rnum,
2649 mcrunh.get_date(),
2650 corsika,
2651 1,
2652 heights,
2653 mcevth_2[0].get_slope(),
2654 mcrunh.get_ELow(),
2655 mcrunh.get_EUpp(),
2656 mcrunh.get_EGS4(),
2657 mcrunh.get_NKG(),
2658 mcrunh.get_Ecutoffh(),
2659 mcrunh.get_Ecutoffm(),
2660 mcrunh.get_Ecutoffe(),
2661 mcrunh.get_Ecutoffg(),
2662 constantC,
2663 constantCKA,
2664 constantCETA,
2665 constantCSTRBA,
2666 constantAATM,
2667 constantBATM,
2668 constantCATM,
2669 constantNFL,
2670 viewcone,
2671 mcrunh.get_wobble(),
2672 mcrunh.get_atmophere()
2673 );
2674
2675 // Store qe for each PMT in output file
2676 TArrayF qe_pmt;
2677 TArrayF wav_pmt;
2678
2679 for(int ict=0;ict<ct_Number;ict++){
2680 McConfigRunHeader[ict]->InitSizePMTs(ct_NPixels);
2681 for(int i=0; i<(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();i++){
2682 McConfigRunHeader[ict]->AddPMT(i);
2683 MGeomPMT &pmt = McConfigRunHeader[ict]->GetPMT(i);
2684 qe_pmt.Set(pointsQE[ict],QE[ict][i][1]);
2685 wav_pmt.Set(pointsQE[ict],QElambda);
2686 pmt.SetArraySize(pointsQE[ict]);
2687 pmt.SetPMTContent(i,wav_pmt,qe_pmt);
2688 }
2689
2690 // Store Light Guides factors in the output file
2691 TArrayF theta_lg;
2692 TArrayF factor_lg;
2693 theta_lg.Set(pointsWC,WC[0]);
2694 factor_lg.Set(pointsWC,WC[1]);
2695 McConfigRunHeader[ict]->SetLightGuides(theta_lg,factor_lg);
2696
2697 }
2698
2699 // Fill the Header Tree with the current leaves of each branch
2700 HeaderTree.Fill() ;
2701
2702 //++
2703 // put the Event to the root file
2704 //--
2705
2706 outfile_temp.Write() ;
2707 outfile_temp.Close() ;
2708
2709 // close input file
2710
2711 if (Trigger_Loop){
2712 log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
2713 ntshow, ntcph[0] );
2714 datafile<<ntshow<<" event(s), with a total of "<<ntcph[0]<<" C.photons"<<endl;
2715 log( SIGNATURE, "Trigger Mode. \n");
2716 log( SIGNATURE, "Fraction of triggers: \n");
2717 datafile<<"Fraction of triggers: "<<endl;
2718 for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
2719 for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
2720 for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
2721 log( SIGNATURE, "Thres %5.1f, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n",
2722 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);
2723 datafile<<"Thres "<<fthrescount<<", Multi "<<imulticount<<", Topo"<<isorttopo[itopocount]<<": ";
2724 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;
2725 }
2726 }
2727 }
2728 }
2729 else{
2730 for(int ict=0;ict<ct_Number;ict++){
2731 log( SIGNATURE,
2732 "%d event(s), with a total of %d C.photons in CT %i (%s)\n",
2733 ntshow, ntcph[ict],ict,GeometryName[ict] );
2734 datafile<<ntshow<<" event(s), with a total of "<<ntcph[ict]
2735 <<" C.photons in CT "<<ict<<" ("<<GeometryName[ict]<<")"<<endl;
2736 log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
2737 ((float)ntrigger[ict]) / ((float)ntshow) * 100.0, ntrigger[ict], ntshow);
2738 datafile<<"Fraction of triggers: "<<((float)ntrigger[ict]) / ((float)ntshow) * 100.0<<" ("<<ntrigger[ict]<<" out of "<<ntshow<<" )"<<endl;
2739 }
2740 }
2741
2742 // close files
2743
2744 log( SIGNATURE, "Closing files\n" );
2745
2746 if( ! Data_From_STDIN ){
2747 for(int ict=0;ict<ct_Number;ict++)
2748 fclose( inputfile[ict] );
2749 }
2750 datafile.close();
2751
2752 // program finished
2753
2754 log( SIGNATURE, "Done.\n");
2755
2756 return( 0 );
2757}
2758//!@}
2759
2760// @T \newpage
2761
2762//!@subsection Functions definition.
2763
2764//!-----------------------------------------------------------
2765// @name present
2766//
2767// @desc Make some presentation
2768//
2769// @date Sat Jun 27 05:58:56 MET DST 1998
2770//------------------------------------------------------------
2771// @function
2772
2773//!@{
2774void
2775present(void)
2776{
2777 cout << "##################################################\n"
2778 << SIGNATURE << '\n' << '\n'
2779 << "Processor of the reflector output\n"
2780 << "J C Gonzalez, Jun 1998\n"
2781 << "##################################################\n\n"
2782 << flush ;
2783}
2784//!@}
2785
2786
2787//!-----------------------------------------------------------
2788// @name usage
2789//
2790// @desc show help
2791//
2792// @date Tue Dec 15 16:23:30 MET 1998
2793//------------------------------------------------------------
2794// @function
2795
2796//!@{
2797void
2798usage(void)
2799{
2800 present();
2801 cout << "\nusage ::\n\n"
2802 << "\t camera "
2803 << " [ -@ paramfile ] "
2804 << " [ -h ] "
2805 << "\n\n or \n\n"
2806 << "\t camera < paramfile"
2807 << "\n\n";
2808 exit(0);
2809}
2810//!@}
2811
2812
2813//!-----------------------------------------------------------
2814// @name log
2815//
2816// @desc function to send log information
2817//
2818// @var funct Name of the caller function
2819// @var fmt Format to be used (message)
2820// @var ... Other information to be shown
2821//
2822// @date Sat Jun 27 05:58:56 MET DST 1998
2823//------------------------------------------------------------
2824// @function
2825
2826//!@{
2827void
2828log(const char *funct, char *fmt, ...)
2829{
2830 va_list args;
2831
2832 // Display the name of the function that called error
2833 printf("[%s]: ", funct);
2834
2835 // Display the remainder of the message
2836 va_start(args, fmt);
2837 vprintf(fmt, args);
2838 va_end(args);
2839}
2840//!@}
2841
2842
2843//!-----------------------------------------------------------
2844// @name error
2845//
2846// @desc function to send an error message, and abort the program
2847//
2848// @var funct Name of the caller function
2849// @var fmt Format to be used (message)
2850// @var ... Other information to be shown
2851//
2852// @date Sat Jun 27 05:58:56 MET DST 1998
2853//------------------------------------------------------------
2854// @function
2855
2856//!@{
2857void
2858error(const char *funct, char *fmt, ...)
2859{
2860 va_list args;
2861
2862 // Display the name of the function that called error
2863 fprintf(stdout, "ERROR in %s: ", funct);
2864
2865 // Display the remainder of the message
2866 va_start(args, fmt);
2867 vfprintf(stdout, fmt, args);
2868 va_end(args);
2869
2870 perror(funct);
2871
2872 exit(1);
2873}
2874//!@}
2875
2876
2877//!-----------------------------------------------------------
2878// @name isA
2879//
2880// @desc returns TRUE(FALSE), if the flag is(is not) the given
2881//
2882// @var s1 String to be searched
2883// @var flag Flag to compare with string s1
2884// @return TRUE: both strings match; FALSE: oth.
2885//
2886// @date Wed Jul 8 15:25:39 MET DST 1998
2887//------------------------------------------------------------
2888// @function
2889
2890//!@{
2891int
2892isA( char * s1, const char * flag ) {
2893 return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
2894}
2895//!@}
2896
2897
2898//!-----------------------------------------------------------
2899// @name read_QE
2900//
2901// @desc read QE data
2902//
2903// @date thu 5 17:59:57 CEST 2002
2904//------------------------------------------------------------
2905// @function
2906
2907//!@{
2908void
2909read_QE(char fname[256], int ict){
2910 ifstream qefile;
2911 char line[LINE_MAX_LENGTH];
2912 int i, j, icount;
2913 float qe;
2914
2915 //------------------------------------------------------------
2916 // second, pixels' QE
2917
2918 // try to open the file
2919
2920 log("read_QE", "Opening the file \"%s\" . . .\n", fname);
2921
2922 qefile.open( fname );
2923
2924 // if it is wrong or does not exist, exit
2925
2926 if ( qefile.bad() )
2927 error( "read_QE", "Cannot open \"%s\". Exiting.\n", fname );
2928
2929 // read file
2930
2931 log("read_QE", "Reading data . . .\n");
2932
2933 i=-1;
2934 icount = 0;
2935
2936 while ( ! qefile.eof() ) {
2937
2938 // get line from the file
2939
2940 qefile.getline(line, LINE_MAX_LENGTH);
2941
2942 // skip if comment
2943
2944 if ( *line == '#' )
2945 continue;
2946
2947 // if it is the first valid value, it is the number of QE data points
2948
2949 if ( i < 0 ) {
2950
2951 // get the number of datapoints
2952
2953 sscanf(line, "%d", &pointsQE[ict]);
2954
2955 // allocate memory for the table of QEs
2956
2957 QE[ict] = new float ** [ct_NPixels];
2958
2959 for ( i=0; i<ct_NPixels; ++i ) {
2960 QE[ict][i] = new float * [2];
2961 QE[ict][i][0] = new float[pointsQE[ict]];
2962 QE[ict][i][1] = new float[pointsQE[ict]];
2963 }
2964
2965 QElambda = new float [pointsQE[ict]];
2966
2967 for ( i=0; i<pointsQE[ict]; ++i ) {
2968 qefile.getline(line, LINE_MAX_LENGTH);
2969 sscanf(line, "%f", &QElambda[i]);
2970 }
2971
2972 i=0;
2973
2974 continue;
2975 }
2976
2977 // get the values (num-pixel, num-datapoint, QE-value)
2978
2979 if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
2980 break;
2981
2982 if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2983 ((j-1) < pointsQE[ict]) && ((j-1) > -1) ) {
2984 QE[ict][i-1][0][j-1] = QElambda[j-1];
2985 QE[ict][i-1][1][j-1] = qe;
2986 }
2987
2988 if ( i > ct_NPixels) break;
2989
2990 icount++;
2991
2992 }
2993
2994 if(icount/pointsQE[ict] < ct_NPixels){
2995 error( "read_QE", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
2996 icount/pointsQE[ict], ct_NPixels );
2997 }
2998
2999 // close file
3000
3001 qefile.close();
3002
3003 // test QE
3004
3005 for(icount=0; icount< ct_NPixels; icount++){
3006 for(i=0; i<pointsQE[ict]; i++){
3007 if( QE[ict][icount][0][i] < 100. || QE[ict][icount][0][i] > 1000. ||
3008 QE[ict][icount][1][i] < 0. || QE[ict][icount][1][i] > 100.){
3009 error( "read_QE", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
3010 icount, i, QE[ict][icount][0][i], QE[ict][icount][1][i] );
3011 }
3012 }
3013 }
3014
3015 // end
3016
3017 log("read_QE", "Done.\n");
3018}
3019//!@}
3020
3021//!-----------------------------------------------------------
3022// @name read_WC
3023//
3024// @desc read WC data
3025//
3026// @date thu 5 17:59:57 CEST 2002
3027//------------------------------------------------------------
3028// @function
3029
3030//!@{
3031void
3032read_WC(void){
3033 ifstream wcfile;
3034 char line[LINE_MAX_LENGTH];
3035 int i;
3036
3037 //------------------------------------------------------------
3038 // Read Light Guides data
3039
3040 // try to open the file
3041
3042 log("read_WC", "Opening the file \"%s\" . . .\n", WC_FILE);
3043
3044 wcfile.open( WC_FILE );
3045
3046 // if it is wrong or does not exist, exit
3047
3048 if ( wcfile.bad() )
3049 error( "read_WC", "Cannot open \"%s\". Exiting.\n", WC_FILE );
3050
3051 // read file
3052
3053 log("read_WC", "Reading data . . .\n");
3054
3055 // get line from the file
3056
3057 wcfile.getline(line, LINE_MAX_LENGTH);
3058
3059 // get the number of datapoints
3060
3061 sscanf(line, "%d", &pointsWC);
3062
3063 // allocate memory for the table of QEs
3064
3065 WC = new float * [2];
3066 WC[0] = new float[pointsWC];
3067 WC[1] = new float[pointsWC];
3068
3069 for ( i=0; i<pointsWC; ++i ) {
3070 wcfile.getline(line, LINE_MAX_LENGTH);
3071 sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
3072 }
3073
3074 // close file
3075
3076 wcfile.close();
3077
3078 // read
3079
3080 log("read_WC", "Done.\n");
3081}
3082//!@}
3083
3084
3085//!-----------------------------------------------------------
3086// @name read_ascii
3087//
3088// @desc read ascii configuration files used by the reflector
3089//
3090// @date tue dec 10 17:14:10 CET 2002
3091//------------------------------------------------------------
3092// @function
3093
3094//!@{
3095void
3096read_ascii(FILE *sp, MMcConfigRunHeader *config)
3097{
3098 Float_t radius = -1.0;
3099 Float_t focal = -1.0;
3100 Float_t point = -1.0;
3101 Float_t spot = -1.0;
3102 Float_t camwidth = -1.0;
3103
3104 Int_t imir;
3105 Float_t f,sx,sy,x,y,z,thetan,phin,xn,yn,zn;
3106 Float_t dx,dy;
3107 Int_t nummir, numref;
3108 Float_t wav,ref;
3109
3110 Char_t token[40];
3111 Char_t line[511];
3112 Char_t flag;
3113
3114 while(1){
3115 if((flag=fgetc(sp))==EOF)
3116 break;
3117 if (flag == '\n') // skip empty lines
3118 continue;
3119
3120 fgets(&line[1],500,sp);
3121 line[0]=flag;
3122
3123 if ( strstr(line, "# reflectivity file") == line ) {
3124 while (strstr(line, "# number of datapoints") != line)
3125 fgets(line,500,sp);
3126
3127 fgets(line,500,sp);
3128 sscanf(line,"%i",&numref);
3129
3130 TArrayF wavarray(numref);
3131 TArrayF refarray(numref);
3132
3133 while (strstr(line, "# datapoints") != line)
3134 fgets(line,500,sp);
3135
3136 for(int i=0; i<numref;i++){
3137 fgets(line,500,sp);
3138 if (line[0] == '#')
3139 {
3140 i--;
3141 continue;
3142 }
3143 sscanf(line,"%f %f",&wav,&ref);
3144 wavarray[i]=wav;
3145 refarray[i]=ref;
3146 }
3147 for (int j=0; j<nummir;j++){
3148
3149 MGeomMirror &mirror = config->GetMirror(j);
3150 mirror.SetArraySize(numref);
3151 mirror.SetReflectivity(wavarray, refarray);
3152 }
3153 continue;
3154 }
3155 if (line[0]== '#')
3156 continue;
3157 if (strstr(line, "type") == line)
3158 continue;
3159 if (strstr(line, "focal_distance") == line){
3160 sscanf(line,"%s %f",token,&focal);
3161 continue;
3162 }
3163 if (strstr(line, "point_spread") == line){
3164 sscanf(line,"%s %f",token,&point);
3165 continue;
3166 }
3167 if (strstr(line, "black_spot") == line){
3168 sscanf(line,"%s %f",token,&spot);
3169 continue;
3170 }
3171 if (strstr(line, "camera_width") == line){
3172 sscanf(line,"%s %f",token,&camwidth);
3173 continue;
3174 }
3175 //
3176 // Skip obsolete magic.def entries:
3177 //
3178 if (strstr(line, "n_pixels") == line)
3179 continue;
3180 if (strstr(line, "pixel_width") == line)
3181 continue;
3182 if (strstr(line, "n_centralpixels") == line)
3183 continue;
3184 if (strstr(line, "n_gappixels") == line)
3185 continue;
3186
3187 if (strstr(line, "n_mirrors") == line){
3188 sscanf(line,"%s %i",token,&nummir);
3189 config->InitSizeMirror(nummir);
3190 continue;
3191 }
3192 if (strstr(line, "r_mirror") == line){
3193 sscanf(line,"%s %f",token,&radius);
3194 continue;
3195 }
3196 if (strstr(line, "define_mirrors") == line){
3197 for(int i=0;i<nummir;i++){
3198 fgets(line,500,sp);
3199 sscanf(line,"%i %f %f %f %f %f %f %f %f %f %f %f",
3200 &imir,&f,&sx,&sy,&x,&y,&z,&thetan,&phin,&xn,&yn,&zn);
3201 config->AddMirror(i);
3202 MGeomMirror &mirror = config->GetMirror(i);
3203 mirror.SetMirrorContent(imir,f,sx,sy,x,y,z,thetan,phin,xn,yn,zn);
3204 }
3205 fgets(line,500,sp);
3206
3207 while ( ! strstr(line, "axis deviation"))
3208 fgets(line,500,sp);
3209
3210 for(int i=0;i<nummir;i++){
3211 fgets(line,500,sp);
3212 if (line[0] == '#')
3213 {
3214 i--;
3215 continue;
3216 }
3217 sscanf(line,"%f %f",&dx,&dy);
3218 MGeomMirror &mirror = config->GetMirror(i);
3219 mirror.SetMirrorDeviations(dx,dy);
3220 }
3221 continue;
3222 }
3223 }
3224 config->SetMagicDef(radius, focal, point, spot, camwidth);
3225}
3226
3227
3228//!-----------------------------------------------------------
3229// @name igen_pixel_coordinates
3230//
3231// @desc generate the pixel center coordinates
3232//
3233// @var *pcam structure camera containing all the
3234// camera information
3235// @return total number of pixels
3236//
3237// DP
3238//
3239// @date Thu Oct 14 10:41:03 CEST 1999
3240//------------------------------------------------------------
3241// @function
3242
3243//!@{
3244/******** igen_pixel_coordinates() *********************************/
3245
3246int igen_pixel_coordinates(struct camera *pcam) {
3247 /* generate pixel coordinates, return value is number of pixels */
3248
3249 int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
3250 float fsegment_fract;
3251 double dtsize;
3252 double dhsize;
3253 double dpsize;
3254 double dxfirst_pix;
3255 double dyfirst_pix;
3256 double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5;
3257 double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5;
3258
3259
3260 double dstartx, dstarty; /* for the gap pixels and outer pixels */
3261 int j, nrow;
3262
3263 dpsize = pcam->dpixdiameter_cm;
3264 dtsize = dpsize * sqrt(3.) / 2.;
3265 dhsize = dpsize / 2.;
3266
3267 /* Loop over central pixels to generate co-ordinates */
3268
3269 for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
3270
3271 /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
3272
3273 pcam->dpixsizefactor[ipixno-1] = 1.;
3274
3275 in = 0;
3276
3277 i = 0;
3278 itot_inside_ring = 0;
3279 iring_no = 0;
3280
3281 /* Calculate the number of pixels out to and including the ring containing pixel number */
3282 /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
3283
3284 while (itot_inside_ring == 0){
3285
3286 iN = 3*(i*(i+1)) + 1;
3287
3288 if (ipixno <= iN){
3289 iring_no = i;
3290 itot_inside_ring = iN;
3291 }
3292
3293 i++;
3294 }
3295
3296
3297 /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
3298
3299 ipix_in_ring = 0;
3300 for (i = 0; i < iring_no; ++i){
3301
3302 ipix_in_ring = ipix_in_ring + 6;
3303 }
3304
3305 /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
3306 /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
3307 /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
3308
3309 isegment = 0;
3310 fsegment_fract = 0.;
3311 if (iring_no > 0) {
3312
3313 isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
3314
3315 fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
3316
3317 }
3318
3319 /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
3320 /* pixel width (flat to flat) dpsize. */
3321
3322 dxfirst_pix = dpsize*iring_no;
3323 dyfirst_pix = 0.;
3324
3325 /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
3326
3327 ddxseg1 = - dhsize*iring_no;
3328 ddyseg1 = dtsize*iring_no;
3329 ddxseg2 = -dpsize*iring_no;
3330 ddyseg2 = 0.;
3331 ddxseg3 = ddxseg1;
3332 ddyseg3 = -ddyseg1;
3333 ddxseg4 = -ddxseg1;
3334 ddyseg4 = -ddyseg1;
3335 ddxseg5 = -ddxseg2;
3336 ddyseg5 = 0.;
3337
3338 /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
3339 /* anti-clockwise around the ring by adding the segment to segment vectors. */
3340
3341 switch (isegment) {
3342
3343 case 0:
3344
3345 pcam->dxc[ipixno-1] = 0.;
3346 pcam->dyc[ipixno-1] = 0.;
3347
3348 case 1:
3349 pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
3350 pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
3351
3352 break;
3353
3354 case 2:
3355
3356 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
3357 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
3358
3359 break;
3360
3361 case 3:
3362
3363 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
3364 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
3365
3366 break;
3367
3368 case 4:
3369
3370 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
3371 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
3372
3373 break;
3374
3375 case 5:
3376
3377 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
3378 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
3379
3380 break;
3381
3382 case 6:
3383
3384 pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
3385 pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
3386
3387 break;
3388
3389 default:
3390
3391 fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
3392 return(0);
3393
3394 } /* end switch */
3395
3396 } /* end for */
3397
3398 dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
3399 dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
3400
3401 if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
3402
3403 j = pcam->inumcentralpixels;
3404
3405 for(i=0; i<pcam->inumgappixels; i=i+6){
3406 pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
3407 pcam->dyc[j + i ] = dstarty;
3408 pcam->dpixsizefactor[j + i] = 1.;
3409 pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
3410 pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
3411 pcam->dpixsizefactor[j + i + 1] = 1.;
3412 pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
3413 pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
3414 pcam->dpixsizefactor[j + i+ 2] = 1.;
3415 pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
3416 pcam->dyc[j + i + 3] = dstarty;
3417 pcam->dpixsizefactor[j + i+ 3] = 1.;
3418 pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
3419 pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
3420 pcam->dpixsizefactor[j + i+ 4] = 1.;
3421 pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
3422 pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
3423 pcam->dpixsizefactor[j + i + 5] = 1.;
3424 } /* end for */
3425 } /* end if */
3426
3427 /* generate positions of the outer pixels */
3428
3429 if( pcam->inumbigpixels > 0 ){
3430
3431 j = pcam->inumcentralpixels + pcam->inumgappixels;
3432
3433 for(i=0; i<pcam->inumbigpixels; i++){
3434 pcam->dpixsizefactor[j + i] = 2.;
3435 }
3436
3437 in = 0;
3438
3439 nrow = (int) ceil(dstartx / 2. / dpsize);
3440
3441 while(in < pcam->inumbigpixels){
3442
3443 pcam->dxc[j + in] = dstartx + dpsize;
3444 pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
3445 pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
3446 pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
3447 pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
3448 pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
3449 pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
3450 pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
3451 pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
3452 pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
3453 pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
3454 pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
3455 for(i=1; i<nrow; i++){
3456 pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
3457 pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
3458 pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
3459 pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
3460 pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
3461 pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
3462 pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
3463 pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
3464 pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
3465 pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
3466 pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
3467 pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
3468 }
3469 in = in + 6 * nrow;
3470 dstartx = dstartx + 2. * dpsize;
3471 nrow = nrow + 1;
3472 } /* end while */
3473
3474 } /* end if */
3475
3476 return(pcam->inumpixels);
3477
3478}
3479//!@}
3480
3481//!-----------------------------------------------------------
3482// @name bpoint_is_in_pix
3483//
3484// @desc check if a point (x,y) in camera coordinates is inside a given pixel
3485//
3486// @var *pcam structure camera containing all the
3487// camera information
3488// @var dx, dy point coordinates in centimeters
3489// @var ipixnum pixel number (starting at 0)
3490// @return TRUE if the point is inside the pixel, FALSE otherwise
3491//
3492// DP
3493//
3494// @date Thu Oct 14 16:59:04 CEST 1999
3495//------------------------------------------------------------
3496// @function
3497
3498//!@{
3499
3500/******** bpoint_is_in_pix() ***************************************/
3501
3502#define sqrt13 0.577350269 // = 1./sqrt(3.)
3503#define sqrt3 1.732050807 // = sqrt(3.)
3504
3505int bpoint_is_in_pix(double dx, double dy, MGeomCam *pgeomcam)
3506{
3507 /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
3508 /* the pixel is assumed to be a "closed set" */
3509
3510 /*
3511 a = length of one of the edges of one pixel,
3512 b = half the width of one pixel
3513 */
3514
3515 const int numN = pgeomcam->GetNumPixels();
3516
3517 for (int i=0; i<numN; i++)
3518 {
3519 MGeomPix &pixel = (*pgeomcam)[i];
3520 const double b = pixel.GetD()/2;
3521 const double a = pixel.GetD()/sqrt3;
3522
3523 const double xx = dx - pixel.GetX();
3524 const double yy = dy - pixel.GetY();
3525
3526 if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
3527 ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){
3528
3529 return i; // inside i
3530 }
3531
3532 // return -1; // outside
3533 }
3534
3535 return -1; // outside
3536}
3537
3538//!@}
3539
3540//------------------------------------------------------------
3541// @name dist_r_P
3542//
3543// @desc distance straight line r - point P
3544//
3545// @date Sat Jun 27 05:58:56 MET DST 1998
3546// @function @code
3547//------------------------------------------------------------
3548// dist_r_P
3549//
3550// distance straight line r - point P
3551//------------------------------------------------------------
3552
3553float
3554dist_r_P(float a, float b, float c,
3555 float l, float m, float n,
3556 float x, float y, float z)
3557{
3558 return (
3559 sqrt((SQR((a-x)*m-(b-y)*l) +
3560 SQR((b-y)*n-(c-z)*m) +
3561 SQR((c-z)*l-(a-x)*n))/
3562 (SQR(l)+SQR(m)+SQR(n))
3563 )
3564 );
3565}
3566
3567//------------------------------------------------------------
3568// @name check_reflector_file
3569//
3570// @desc check if a given reflector file has the right signature
3571// @desc return TRUE or FALSE
3572//
3573// @date Mon Feb 14 16:44:21 CET 2000
3574// @function @code
3575//------------------------------------------------------------
3576
3577int check_reflector_file(FILE *infile){
3578
3579 char sign[20]; // auxiliary variable
3580
3581 strcpy(sign, REFL_SIGNATURE_B);
3582
3583 fread( (char *)sign, strlen(REFL_SIGNATURE_B), 1, infile);
3584 if (strcmp(sign, REFL_SIGNATURE_A) == 0){
3585 fread( (char *)sign, 1, 1, infile);
3586 return 4;
3587 }
3588 else if (strcmp(sign, REFL_SIGNATURE_B) == 0){
3589 fread( (char *)sign, 1, 1, infile);
3590 return 5;
3591 }
3592 else if (strcmp(sign, REFL_SIGNATURE_C) == 0){
3593 // An empty bin has been removed and therefore we do not need to rezd it.
3594 return 6;
3595 }
3596 else {
3597 cout << "ERROR: Signature of .rfl file is not correct\n";
3598 cout << '"' << sign << '"' << '\n';
3599 cout << "should be: " << REFL_SIGNATURE_A <<" or "<< REFL_SIGNATURE_B <<" or " << REFL_SIGNATURE_C <<" or "<< '\n';
3600 return(FALSE);
3601 }
3602
3603
3604}
3605
3606//------------------------------------------------------------
3607// @name lin_interpol
3608//
3609// @desc interpolate linearly between two points returning the
3610// @desc y-value of the result
3611//
3612// @date Thu Feb 17 11:31:32 CET 2000
3613// @function @code
3614//------------------------------------------------------------
3615
3616float lin_interpol( float x1, float y1, float x2, float y2, float x){
3617
3618 if( (x2 - x1)<=0. ){ // avoid division by zero, return average
3619 cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
3620 return((y1+y2)/2.);
3621 }
3622 else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
3623 return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
3624 }
3625}
3626
3627//------------------------------------------------------------
3628// @name size_rotated
3629//
3630// @desc It rotates the NSB
3631//
3632// @date Tue Jan 7 17:09:25 CET 2003
3633// @function @code
3634//------------------------------------------------------------
3635
3636int size_rotated(
3637 float *rotated,
3638 float nsb[iMAXNUMPIX],
3639 float rho)
3640{
3641 int r=0;
3642 float size_rotated;
3643 float Num_Pixels;
3644 float PixNum=0;
3645 float rho_pixel;
3646 int j=0;
3647 int k=0;
3648
3649 for(int i=1; i<iMAXNUMPIX;i++){
3650// Substituir per Int_t radius[iMAXNUMPIX]={0,1,1,1,1,1,1,2,2,2,2, ...}
3651 Num_Pixels=0;
3652 for (int ii=1; ii<17;ii++){
3653 if (Num_Pixels >= i){
3654 r=ii-1;
3655 break;
3656 }
3657
3658
3659
3660 if (ii<12){
3661 Num_Pixels+=ii*6;
3662 PixNum=6*ii;
3663 // size_rotated = (360/PixNum);
3664 //rho_pixel = rho/size_rotated;
3665 }
3666
3667
3668 else
3669 {
3670 Num_Pixels+=(ii-6)*6;
3671 PixNum=(ii-6)*6;
3672 // size_rotated = (360/PixNum);
3673 // rho_pixel = rho/size_rotated;
3674 }
3675 }
3676
3677
3678
3679 size_rotated = (360/PixNum);
3680 rho_pixel = rho/size_rotated;
3681
3682 //Buscar j i k que no canvin de r!!!
3683
3684 j=i-int(rho_pixel);
3685
3686 //cout<<"j inicial "<<j<<endl;
3687 int MINPIXinner[13]= {0,1,7,19,37,61,91,127,169,217,271,331,397};
3688 int MINPIXouter[5]={397,433,475,523,578};
3689
3690 if( r<12 ){
3691 if (j < MINPIXinner[r])
3692 {
3693 j = i+6 * r - int(rho_pixel);
3694 }
3695
3696 k=j-1;
3697 if (k < MINPIXinner[r])
3698 {
3699 k = MINPIXinner[r+1]-1;
3700 }
3701 }
3702 if( r > 11 && r < 16){
3703 if (j < MINPIXouter[r-12])
3704 {
3705 j = i + (r - 6) * 6 - int(rho_pixel);
3706
3707 }
3708
3709 k=j-1;
3710 if ( k < MINPIXouter[r-12])
3711 {
3712 k = MINPIXouter[r-11] - 1;
3713 }
3714 }
3715 rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k];
3716 }
3717 return (int(rho_pixel));
3718
3719}
3720
3721
3722//------------------------------------------------------------
3723// @name produce_phes
3724//
3725// @desc read the photons of an event, pixelize them and simulate QE
3726// @desc return various statistics and the array of Photoelectrons
3727//
3728// @date Mon Feb 14 16:44:21 CET 2000
3729// @function @code
3730//------------------------------------------------------------
3731
3732int produce_phes( FILE *sp, // the input file
3733 class MGeomCam *camgeom, // the camera layout
3734 float minwl_nm, // the minimum accepted wavelength
3735 float maxwl_nm, // the maximum accepted wavelength
3736 class MTrigger *trigger, // the generated phes
3737 class MFadc *fadc,
3738 int *itotnphe, // total number of produced photoelectrons
3739 float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
3740 int *incph, // total number of cph within camera
3741 float *tmin_ns, // minimum arrival time of all phes
3742 float *tmax_ns, // maximum arrival time of all phes
3743 int ict // Telescope that is being analised to get the right QE.
3744 ){
3745
3746 static uint i;
3747 static int k, ipixnum;
3748 static float cx, cy, wl, qe, t;
3749 static float cu, cv, cw;
3750 static MCCphoton photon;
3751 static float **qept;
3752 static char flag[SIZE_OF_FLAGS + 1];
3753 static float radius_mm;
3754
3755 // reset variables
3756
3757 for ( i=0; i<camgeom->GetNumPixels(); ++i ){
3758
3759 nphe[i] = 0.0;
3760
3761 }
3762
3763 *incph = 0;
3764
3765 radius_mm = camgeom->GetMaxRadius();
3766
3767 TRandom random;
3768 random.SetSeed((Int_t)(get_seeds(0)*get_seeds(1)));
3769
3770 float spotsigma = get_sigma_xy_cm_spot();
3771
3772 //- - - - - - - - - - - - - - - - - - - - - - - - -
3773 // read photons and "map" them into the pixels
3774 //--------------------------------------------------
3775
3776 // initialize CPhoton
3777
3778 photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
3779
3780 // read the photons data
3781
3782
3783 // loop over the photons
3784
3785 while (1) {
3786
3787 fread ( flag, SIZE_OF_FLAGS, 1, sp );
3788
3789 if (isA( flag, FLAG_END_OF_EVENT ))
3790 break;
3791
3792 memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
3793
3794 fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
3795
3796
3797 // Check if photon is inside trigger time range
3798
3799 t = photon.get_t() ;
3800
3801 if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
3802 continue;
3803 }
3804
3805 //
3806 // Pixelization
3807 //
3808
3809 cx = photon.get_x();
3810 cy = photon.get_y();
3811
3812
3813 if (spotsigma > 0.)
3814 {
3815 cx += random.Gaus(0.,spotsigma);
3816 cy += random.Gaus(0.,spotsigma);
3817 }
3818
3819 // get wavelength
3820
3821 wl = photon.get_wl();
3822
3823 // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n";
3824
3825 // check if photon has valid wavelength and is inside outermost camera radius
3826
3827 if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy)*10 > radius_mm ) ){
3828 continue;
3829
3830 }
3831
3832 ipixnum = bpoint_is_in_pix(cx*10, cy*10, camgeom);
3833
3834 // -1 = the photon is in none of the pixels
3835 // 0 = the phton is in the central pixel, which is not used for trigger
3836 if (ipixnum==-1 || ipixnum==0) {
3837 continue;
3838 }
3839
3840 // AM changed meaning of incph: before it was all photons read from
3841 // reflector file, now only those within a valid pixel:
3842 //
3843 // increase number of photons
3844 (*incph)++;
3845
3846 //+++
3847 // QE simulation
3848 //---
3849
3850 // set pointer to the QE table of the relevant pixel
3851
3852 qept = (float **)QE[ict][ipixnum];
3853
3854 // check if wl is inside table; outside the table, QE is assumed to be zero
3855
3856 if((wl < qept[0][0]) || (wl > qept[0][pointsQE[ict]-1])){
3857 continue;
3858
3859 }
3860
3861 // find data point in the QE table (-> k)
3862
3863 k = 1; // start at 1 because the condition was already tested for 0
3864 while (k < pointsQE[ict]-1 && qept[0][k] < wl){
3865 k++;
3866 }
3867
3868 // calculate the qe value between 0. and 1.
3869
3870 qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
3871
3872 // Apply incient angular correction due to Light Guides
3873
3874 cu=photon.get_u();
3875 cv=photon.get_v();
3876 cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*3.14159/180.0;
3877
3878 // find data point in the QE table (-> k)
3879
3880 k = 0;
3881 while (k < pointsWC-1 && WC[0][k] < cw){
3882 k++;
3883 }
3884
3885 // correct the qe with WC data.
3886
3887 qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
3888
3889 // if random > quantum efficiency, reject it
3890
3891 if ( (RandomNumber) > qe ) {
3892 continue;
3893 }
3894
3895 //+++
3896 // The photon has produced a photo electron
3897 //---
3898
3899 // cout << " accepted\n";
3900
3901 // increment the number of photoelectrons in the relevant pixel
3902
3903 nphe[ipixnum] += 1.0;
3904
3905 // store the new photoelectron
3906
3907 fadc->Fill(ipixnum,(t-*tmin_ns) ,
3908 trigger->FillShow(ipixnum,t-*tmin_ns),
3909 !((*camgeom)[ipixnum].GetD()>(*camgeom)[0].GetD()));
3910
3911 *itotnphe += 1;
3912 }
3913
3914 return(0);
3915
3916}
3917
3918
3919//------------------------------------------------------------
3920// @name produce_nsbrates
3921//
3922// @desc read the starfield file, call produce_phes on it in,
3923// @desc different wavebands, calculate the nsbrates
3924//
3925// @date Mon Feb 14 16:44:21 CET 2000
3926// @function @code
3927//------------------------------------------------------------
3928
3929int produce_nsbrates( char *iname, // the starfield input file name
3930 MGeomCam *camgeom, // camera layout
3931 float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS], // the product of this function:
3932 // the NSB rates in phe/ns for each pixel
3933 int ict
3934 ){
3935
3936 uint i, j; // counters
3937 int k, ii; // counters
3938
3939 static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
3940 WAVEBANDBOUND2,
3941 WAVEBANDBOUND3,
3942 WAVEBANDBOUND4,
3943 WAVEBANDBOUND5,
3944 WAVEBANDBOUND6 };
3945
3946 FILE *infile; // the input file
3947 fpos_t fileposition; // marker on the input file
3948 static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
3949 static MCEventHeader evth; // the event header
3950 static MCEventHeader evth_2; // the event header
3951 static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
3952 int reflector_file_version;
3953 int itnphe; // total number of produced photoelectrons
3954 int itotnphe; // total number of produced photoelectrons after averaging
3955 int incph; // total number of cph read
3956 float tmin_ns; // minimum arrival time of all phes
3957 float tmax_ns; // maximum arrival time of all phes
3958 float integtime_ns; // integration time from the starfield generator
3959 char flag_new[4];
3960
3961 // open input file
3962
3963 log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
3964
3965 infile = fopen( iname, "r" );
3966
3967 // check if the satrfield input file exists
3968
3969 if ( infile == NULL ) {
3970
3971 log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
3972
3973 log( SIGNATURE, "There is not NSB from the Stars\n");
3974
3975 return (0);
3976 }
3977
3978 // get signature, and check it
3979
3980 if((reflector_file_version=check_reflector_file( infile ))==FALSE){
3981 exit(1);
3982 }
3983
3984 // Instance of MTrigger and MFadc
3985
3986 MTrigger trigger(camgeom->GetNumPixels(),camgeom,Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
3987 MFadc flashadc(camgeom->GetNumPixels());
3988
3989 // initialize flag
3990
3991 strcpy( cflag, " \0" );
3992
3993 // get flag
3994
3995 fread( cflag, SIZE_OF_FLAGS, 1, infile );
3996
3997 if ( ! feof(infile)){
3998
3999 // reading .rfl file
4000
4001 if(!isA( cflag, FLAG_START_OF_RUN )){
4002 error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
4003 }
4004 else { // found start of run
4005
4006 if (reflector_file_version > 5){
4007
4008 fread( flag_new, 4, 1, infile );
4009
4010 if(!isA( flag_new, FLAG_START_OF_HEADER)){
4011
4012 // We break the main loop
4013 cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
4014 cout<<" We assume no Star Light"<<endl;
4015 return(0);
4016 }
4017
4018 Float_t flag_temp[SIZE_OF_MCRUNHEADER];
4019
4020 fread( flag_temp, (SIZE_OF_MCRUNHEADER)*sizeof(float), 1, infile );
4021
4022 }
4023
4024 fread( cflag, SIZE_OF_FLAGS, 1, infile );
4025
4026 if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
4027 fread( flag_new, 4, 1, infile );
4028
4029 if(!isA( flag_new, FLAG_EVENT_HEADER)){
4030
4031 // We break while events loop
4032 cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<endl;
4033 cout<<" We assume no light from Stars"<<endl;
4034 return(0);
4035 }
4036
4037 // get MCEventHeader
4038
4039 if (reflector_file_version<6)
4040 fread( (char*)&evth, evth.mysize(), 1, infile );
4041 else
4042 fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
4043
4044 if (reflector_file_version<6)
4045 integtime_ns = evth.get_energy();
4046 else
4047 integtime_ns = evth_2.get_energy();
4048
4049 // memorize where we are in the file
4050
4051 if (fgetpos( infile, &fileposition ) != 0){
4052 error( SIGNATURE, "Cannot position in file ...\n");
4053 }
4054
4055 // loop over the wavebands
4056
4057 for(i=0; i<iNUMWAVEBANDS; i++){
4058
4059 // initialize the rate array
4060
4061 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4062 rate_phepns[j][i] = 0.;
4063 }
4064
4065 itotnphe = 0;
4066
4067 // read the photons and produce the photoelectrons
4068 // - in order to average over the QE simulation, call the
4069 // production function iNUMNSBPRODCALLS times
4070
4071 for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
4072
4073 // position the file pointer to the beginning of the photons
4074
4075 fsetpos( infile, &fileposition);
4076
4077 // produce photoelectrons
4078
4079 k = produce_phes( infile,
4080 camgeom,
4081 wl_nm[i],
4082 wl_nm[i+1],
4083 &trigger, // this is a dummy here
4084 &flashadc, // this is a dummy here
4085 &itnphe,
4086 nphe, // we want this!
4087 &incph,
4088 &tmin_ns,
4089 &tmax_ns,
4090 0
4091 ); // photons from starfield
4092
4093 if( k != 0 ){ // non-zero returnvalue means error
4094 cout << "Exiting.\n";
4095 exit(1);
4096 }
4097
4098 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
4099 rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
4100 }
4101
4102 itotnphe += itnphe;
4103
4104 } // end for(ii=0 ...
4105
4106 fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
4107 wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
4108
4109 } // end for(i=0 ...
4110
4111 }
4112 else{
4113 cout << "Starfield file contains no event.\nExiting.\n";
4114 exit(1);
4115 } // end if( isA ... event
4116 } // end if ( isA ... run
4117 }
4118 else{
4119 cout << "Starfield file contains no run.\nExiting.\n";
4120 exit(1);
4121 }
4122
4123 fclose( infile );
4124
4125 return(0);
4126
4127}
4128
4129
4130//------------------------------------------------------------
4131// @name produce_nsb_phes
4132//
4133// @desc produce the photoelectrons from the nsbrates
4134//
4135// @date Thu Feb 17 17:10:40 CET 2000
4136// @function @code
4137//------------------------------------------------------------
4138
4139int produce_nsb_phes( float *atmin_ns,
4140 float *atmax_ns,
4141 float theta_rad,
4142 struct camera *cam,
4143 float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
4144 float difnsbr_phepns[iMAXNUMPIX],
4145 float extinction[iNUMWAVEBANDS],
4146 float fnpx[iMAXNUMPIX],
4147 MTrigger *trigger,
4148 MFadc *fadc,
4149 int *inphe,
4150 float base_mv[iMAXNUMPIX]){
4151
4152 float simtime_ns; // NSB simulation time
4153 int i, j, k, ii;
4154 float zenfactor; // correction factor calculated from the extinction
4155 int inumnsbphe; // number of photoelectrons caused by NSB
4156
4157 float t;
4158 TRandom random; // Random numbers generator
4159
4160 random.SetSeed(get_seeds(1));
4161
4162 ii = *inphe; // avoid dereferencing
4163
4164 // check if the arrival times are set; if not generate them
4165
4166 if(*atmin_ns <SIMTIMEOFFSET_NS || *atmin_ns > *atmax_ns){
4167 *atmin_ns = 0.;
4168 *atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
4169
4170 }
4171 else{ // extend the event time window by the given offsets
4172
4173 *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
4174 *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
4175
4176 simtime_ns = *atmax_ns - *atmin_ns;
4177
4178 // make sure the simulated time is long enough for the FADC
4179 // simulation and not too long
4180
4181 if(simtime_ns< TOTAL_TRIGGER_TIME){
4182 *atmin_ns = *atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
4183 *atmax_ns = *atmin_ns + TOTAL_TRIGGER_TIME;
4184 simtime_ns = TOTAL_TRIGGER_TIME;
4185 }
4186
4187 if(simtime_ns> TOTAL_TRIGGER_TIME){
4188 *atmax_ns =*atmin_ns + TOTAL_TRIGGER_TIME;
4189 simtime_ns = TOTAL_TRIGGER_TIME;
4190 }
4191
4192 }
4193
4194 // initialize baselines
4195
4196 for(i=0; i<cam->inumpixels; i++){
4197 base_mv[i] = 0.;
4198 }
4199
4200 // calculate baselines and generate phes
4201
4202 for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
4203
4204 // calculate the effect of the atmospheric extinction
4205
4206 zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
4207
4208 for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
4209
4210 inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
4211 * simtime_ns );
4212
4213 base_mv[j] += inumnsbphe;
4214
4215 // randomize
4216
4217 if (inumnsbphe>0.0){
4218 inumnsbphe = random.Poisson (inumnsbphe );
4219 }
4220
4221 // create the photoelectrons
4222
4223 for(k=0; k < inumnsbphe; k++){
4224
4225 t=(RandomNumber * simtime_ns);
4226
4227 (*fadc).Fill(j,t ,(*trigger).FillNSB(j,t));
4228
4229 ii++; // increment total number of photoelectons
4230
4231 fnpx[j] += 1.; // increment number of photoelectrons in this pixel
4232
4233 }
4234
4235 } // end for(j=0 ...
4236 } // end for(i=0 ...
4237
4238 // finish baseline calculation
4239
4240 for(i=0; i<cam->inumpixels; i++){
4241 base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
4242 }
4243
4244 *inphe = ii; // update the pointer
4245
4246 return(0);
4247}
4248
4249
4250// @endcode
4251
4252
4253//=------------------------------------------------------------
4254//!@subsection Log of this file.
4255
4256//!@{
4257//
4258// $Log: not supported by cvs2svn $
4259// Revision 1.67 2004/01/30 09:51:18 blanch
4260// [Changes mainly done by A. Moralejo]
4261//
4262// Several minnor changes have been done. For instance, some name of the
4263// variables have been modified to a more self-explained name and
4264// modifications while reading the asciis files at the end of the reflector file
4265// has been introduced.
4266//
4267// The possibilty to unlarge the point spread function has been introduced
4268// in order to be able to simualte the current data.
4269//
4270// All pixels are always written.
4271//
4272// Revision 1.65 2003/10/26 19:43:00 blanch
4273// - The screen output information has been improved to prevent some
4274// non-desired running conditions just looking at the output screen.
4275// - One MMcEvt for each Telscope is stored in the output file.
4276// - 500 empty events are simualted to get a more precise estimation of the
4277// pedestal Sigma for each pixel.
4278//
4279// Revision 1.64 2003/10/21 07:42:50 blanch
4280// A factor 2.35 to transform the fwhm into the sigma of gaussian was missing
4281// in the storing of FADC single hpe pulse determination.
4282//
4283// Revision 1.63 2003/10/17 19:38:31 blanch
4284// Now the camera program will stop if a undefined Geometry is required.
4285// The NSB is internally scaled for any camera geometry and qe.
4286// The seeds in the input card are used to initilize the random numbers.
4287// The Amplitud stored in the MMcFadcHeader is the amplitud of the sphe reponse.
4288// The Pedestal rms is simulated in an artificial empty event.
4289//
4290// Revision 1.62 2003/09/26 11:25:07 blanch
4291// Modification to be able to read MGeomCam branch for any Geometry.
4292//
4293// Revision 1.61 2003/09/25 17:09:20 blanch
4294// Bug on the number of phe from diffuse NSB fixed.
4295//
4296// Revision 1.60 2003/09/23 16:50:55 blanch
4297// WE do not read ct_file anymore since all Telescope information is
4298// in the reflector or in MGeomCam.
4299//
4300// Revision 1.58 2003/09/15 09:59:53 blanch
4301// The concept of the camera prgoram has not changed but this version has
4302// quite a lot of changes to allow several Camera geometries as well as
4303// multitelescope.
4304//
4305// It is suposed to be the first working code for camera 0.7.
4306//
4307// Revision 1.57 2003/07/17 18:02:46 blanch
4308// Several new features introduced as well as fixed bugs
4309//
4310// - 1/100 events printed out
4311// - Low gain implemented
4312// - Different response for outer and inner pixels
4313// - Some warnings removed
4314// - pedestal and qe file from inpuit card
4315// - Faster electronic simulation
4316//
4317// Revision 1.55 2003/02/12 12:22:10 blanch
4318// *** empty log message ***
4319//
4320// Revision 1.54 2003/02/12 11:55:01 blanch
4321// *** empty log message ***
4322//
4323// Revision 1.53 2003/01/23 18:35:21 blanch
4324// *** empty log message ***
4325//
4326// Revision 1.52 2003/01/20 17:19:20 blanch
4327// It fills the MMcCorsikaRun.
4328//
4329// Revision 1.51 2003/01/14 13:40:17 blanch
4330// MRawRunHeader::fNumEvents has been filled with the number of events in
4331// this file.
4332// Problems in fImpact computation have been solved.
4333// Option to set a dc value to rise the discriminator threshold has been added.
4334//
4335// Revision 1.50 2003/01/07 16:33:31 blanch
4336// Star Field Rotation has been implemented by Raul Orduna. Now there is a
4337// rotation for each shower. It is done by a non enter pixel rotation assuming
4338// a circular symetry of the camera. It is not exact but it is accurate enough and
4339// much faster.
4340//
4341// Revision 1.49 2002/12/13 10:04:07 blanch
4342// *** empty log message ***
4343//
4344// Revision 1.48 2002/12/12 17:40:50 blanch
4345// *** empty log message ***
4346//
4347// Revision 1.47 2002/12/10 17:19:31 blanch
4348// *** empty log message ***
4349//
4350// Revision 1.46 2002/11/08 17:51:00 blanch
4351// *** empty log message ***
4352//
4353// Revision 1.45 2002/10/29 17:15:27 blanch
4354// Reading several reflector versions.
4355//
4356// Revision 1.44 2002/10/18 16:53:03 blanch
4357// Modification to read several reflector files.
4358//
4359// Revision 1.43 2002/09/13 10:53:39 blanch
4360// Minor change to remove some undisired comments.
4361//
4362// Revision 1.42 2002/09/09 16:00:49 blanch
4363// Statement has been included to avoid writting to disk MParContainer and MArray.
4364// It has also been added the effect of the WC, the actual values must be added,
4365// once they are measured.
4366//
4367// Revision 1.41 2002/09/04 09:57:42 blanch
4368// Modifications done to use MGeomCam from MARS.
4369//
4370// Revision 1.40 2002/07/16 16:15:22 blanch
4371// A first implementation of the Star field rotation has been done.
4372//
4373// Revision 1.39 2002/04/27 10:48:39 blanch
4374// Some unused varibles have been removed.
4375//
4376// Revision 1.38 2002/03/18 18:44:29 blanch
4377// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
4378//
4379// Revision 1.37 2002/03/18 16:42:20 blanch
4380// The data member fProductionSite of the MMcRunHeader has been set to 0,
4381// which means that the prodution site is unknown.
4382//
4383// Revision 1.35 2002/03/15 15:15:52 blanch
4384// Several modifications to simulate the actual trigger zone.
4385//
4386// Revision 1.34 2002/03/13 18:13:56 blanch
4387// Some changes to fill correctly the new format of MMcRunHeader.
4388//
4389// Revision 1.33 2002/03/04 17:21:48 blanch
4390// Small and not important changes.
4391//
4392// Revision 1.32 2002/02/28 15:04:52 blanch
4393// A small back has been solved. Before, while not using the option
4394// writte_all_images, not all triggered showers were stored. Now it is solved.
4395// For that it is importatn taht the less restrictive trigger option is
4396// checked first.
4397// A new facility has been introduced and now one can choose the step size in
4398// trigger loop mode for the discriminator threshold.
4399// The close-compact topology for less than 3 pixels does not make sense. Before
4400// the program was ignoring that, now it switch to simple neighbour condition.
4401//
4402// Revision 1.31 2002/01/18 17:41:02 blanch
4403// The option of adding noise to all pixels or to not adding the noise
4404// has been added.
4405// We removed the pixels larger than 577. When there were more than one
4406// trigger in one shower, the pixel number was increasing. Now it is
4407// flagged by the variable MMcTrig::fFirstLvlTrig.
4408//
4409// Revision 1.30 2001/11/27 09:49:54 blanch
4410// Fixing bug which was treating wrongly the extension of star photons.
4411//
4412// Revision 1.29 2001/11/14 17:38:23 blanch
4413// Sveral changes have been done:
4414// - bpoint_in_in_pixel has been dodified to speed up the program
4415// - Some minor changes have been done to remove warnings
4416// - buffer size and split version of the Branches have been removed
4417// - Some modifications were needed to adat the program to the new
4418// MRawEvtData::DeletePixels
4419//
4420// Revision 1.28 2001/10/26 16:31:45 blanch
4421// Removing several warnings.
4422//
4423// Revision 1.27 2001/09/05 10:04:33 blanch
4424// *** empty log message ***
4425//
4426// Revision 1.26 2001/07/19 09:29:53 blanch
4427// Different threshold for each pixel can be used.
4428//
4429// Revision 1.25 2001/05/08 08:07:54 blanch
4430// New numbering for branches from different trigger conditions has been
4431// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
4432// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
4433// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
4434// information already filled.
4435//
4436// Revision 1.24 2001/04/06 16:48:09 magicsol
4437// New camera version able to read the new format of the reflector output:
4438// reflector 0.4
4439//
4440// Revision 1.23 2001/03/28 16:13:41 blanch
4441// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
4442// been solved.
4443//
4444// Revision 1.22 2001/03/20 18:52:43 blanch
4445// *** empty log message ***
4446//
4447// Revision 1.21 2001/03/19 19:53:03 blanch
4448// Some print outs have been removed.
4449//
4450// Revision 1.20 2001/03/19 19:30:06 magicsol
4451// Minor changes have been done to improve the FADC pedestals treatment.
4452//
4453// Revision 1.19 2001/03/05 11:14:41 magicsol
4454// I changed the position of readinf a parameter. It is a minnor change.
4455//
4456// Revision 1.18 2001/03/05 10:36:52 blanch
4457// A branch with information about the FADC simulation (MMcFadcHeader) is writen
4458// in the McHeaders tree of the aoutput root file.
4459// The NSB contribution is only applied if the the number of phe form the shower
4460// are above a given number.
4461//
4462// Revision 1.17 2001/02/23 11:05:57 magicsol
4463// Small changes due to slightly different output format and the introduction of
4464// pedesals for teh FADC.
4465//
4466// Revision 1.16 2001/01/18 18:44:40 magicsol
4467// *** empty log message ***
4468//
4469// Revision 1.15 2001/01/17 09:32:27 harald
4470// The changes are neccessary to have the same name for trees in MC and in
4471// data. So now there should be no differences in MC data and real data from
4472// FADC system.
4473//
4474// Revision 1.14 2001/01/15 12:33:34 magicsol
4475// Some modifications have been done to use the new (Dec'2000) Raw data format.
4476// There are also some minnors modifications to adapt the improvements in the
4477// MTrigger class (overlaping time and trigger cells).
4478//
4479// Revision 1.13 2000/10/25 08:14:23 magicsol
4480// The routine that produce poisson random numbers to decide how many phe
4481// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
4482// is used.
4483//
4484// Revision 1.10 2000/07/04 14:10:20 MagicSol
4485// Some changes have been done in the root output file. The RawEvt tree is only
4486// stored in the single trigger mode.
4487// The trigger input parameters are also given by the general input card.
4488// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
4489// each one can be studied seperately.
4490//
4491// Revision 1.9 2000/06/13 13:25:24 blanch
4492// The multiple arrays have been replaced, since they do not work
4493// in alpha machines. Now we are using pointers and the command new
4494// to allocate memory.
4495//
4496// Revision 1.8 2000/05/11 13:57:27 blanch
4497// The option to loop over several trigger configurations has been included.
4498// Some non-sense with triggertime range has been solved.
4499// Montecarlo information and ADC counts are saved in a root file.
4500// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
4501//
4502// Revision 1.7 2000/03/24 18:10:46 blanch
4503// A first FADC simulation and a trigger simulation are already implemented.
4504// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
4505// 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.
4506//
4507// Revision 1.6 2000/03/20 18:35:11 blanch
4508// 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.
4509//
4510// Revision 1.5 2000/02/18 17:40:35 petry
4511// This version includes drastic changes compared to camera.cxx 1.4.
4512// It is not yet finished and not immediately useful because the
4513// trigger simulation is not yet re-implemented. I had to take it
4514// out together with some other stuff in order to tidy the whole
4515// program up. This is not meant as an insult to anyone. I needed
4516// to do this in order to be able to work on it.
4517//
4518// This version has been put in the repository in order to be
4519// able to share the further development with others.
4520//
4521// If you need something working, wait or take an earlier one.
4522// See file README.
4523//
4524// Revision 1.4 2000/01/25 08:36:23 petry
4525// The pixelization in previous versions was buggy.
4526// This is the first version with a correct pixelization.
4527//
4528// Revision 1.3 2000/01/20 18:22:17 petry
4529// Found little bug which makes camera crash if it finds a photon
4530// of invalid wavelength. This bug is now fixed and the range
4531// of valid wavelengths extended to 290 - 800 nm.
4532// This is in preparation for the NSB simulation to come.
4533// Dirk
4534//
4535// Revision 1.2 1999/11/19 08:40:42 harald
4536// Now it is possible to compile the camera programm under osf1.
4537//
4538// Revision 1.1.1.1 1999/11/05 11:59:31 harald
4539// This the starting point for CVS controlled further developments of the
4540// camera program. The program was originally written by Jose Carlos.
4541// But here you can find a "rootified" version to the program. This means
4542// that there is no hbook stuff in it now. Also the output of the
4543// program changed to the MagicRawDataFormat.
4544//
4545// The "rootification" was done by Dirk Petry and Harald Kornmayer.
4546//
4547// Revision 1.3 1999/10/22 15:01:28 petry
4548// version sent to H.K. and N.M. on Fri Oct 22 1999
4549//
4550// Revision 1.2 1999/10/22 09:44:23 petry
4551// first synthesized version which compiles and runs without crashing;
4552//
4553// Revision 1.1.1.1 1999/10/21 16:35:10 petry
4554// first synthesised version
4555//
4556// Revision 1.13 1999/03/15 14:59:05 gonzalez
4557// camera-1_1
4558//
4559// Revision 1.12 1999/03/02 09:56:10 gonzalez
4560// *** empty log message ***
4561//
4562//
4563//!@}
4564
4565//=EOF
Note: See TracBrowser for help on using the repository browser.