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

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