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

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