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.6 $
|
---|
24 | // $Author: blanch $
|
---|
25 | // $Date: 2000-03-20 18:35:11 $
|
---|
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 "TFile.h"
|
---|
52 | #include "TTree.h"
|
---|
53 | #include "TBranch.h"
|
---|
54 | #include "TCanvas.h"
|
---|
55 |
|
---|
56 | #include "MDiag.h"
|
---|
57 |
|
---|
58 | #include "MRawEvt.h"
|
---|
59 | #include "MMcEvt.h"
|
---|
60 | #include "MMcTrig.hxx"
|
---|
61 |
|
---|
62 | #include "MTrigger.hxx"
|
---|
63 |
|
---|
64 | /*!@"
|
---|
65 |
|
---|
66 | All the defines are located in the file |camera.h|.
|
---|
67 |
|
---|
68 | @"*/
|
---|
69 |
|
---|
70 | #include "camera.h"
|
---|
71 | //!@}
|
---|
72 |
|
---|
73 | /*!@"
|
---|
74 |
|
---|
75 | The following set of flags are used in time of compilation. They do
|
---|
76 | not affect directly the behaviour of the program at run-time
|
---|
77 | (though, of course, if you disconnected the option for
|
---|
78 | implementation of the Trigger logic, you will not be able to use any
|
---|
79 | trigger at all. The 'default' values mean default in the sense of
|
---|
80 | what you got from the server when you obtained this program.
|
---|
81 |
|
---|
82 | @"*/
|
---|
83 |
|
---|
84 | //!@{
|
---|
85 |
|
---|
86 | // flag for debugging (default: OFF )
|
---|
87 | #define __DEBUG__
|
---|
88 | #undef __DEBUG__
|
---|
89 |
|
---|
90 | //!@}
|
---|
91 |
|
---|
92 | //=-----------------------------------------------------------
|
---|
93 | //!@subsection Definition of global variables.
|
---|
94 |
|
---|
95 | /*!@"
|
---|
96 |
|
---|
97 | Now we define some global variables with data about the telescope,
|
---|
98 | such as "focal distance", number of pixels/mirrors,
|
---|
99 | "size of the camera", and so on.
|
---|
100 |
|
---|
101 | @"*/
|
---|
102 |
|
---|
103 | /*!@"
|
---|
104 |
|
---|
105 | Depending on the telescope we are using (CT1 or MAGIC), the
|
---|
106 | information stored in the definition file is different.
|
---|
107 | The variable |ct_Type| has the value 0 when we use
|
---|
108 | CT1, and 1 when we use MAGIC.
|
---|
109 |
|
---|
110 | @"*/
|
---|
111 |
|
---|
112 | //!@{
|
---|
113 | static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC
|
---|
114 | //!@}
|
---|
115 |
|
---|
116 | /*!@"
|
---|
117 |
|
---|
118 | And this is the information about the whole telescope.
|
---|
119 |
|
---|
120 | @"*/
|
---|
121 |
|
---|
122 | //!@{
|
---|
123 |
|
---|
124 | // parameters of the CT (from the CT definition file)
|
---|
125 |
|
---|
126 | ////@: Focal distances [cm]
|
---|
127 | //static float *ct_Focal;
|
---|
128 |
|
---|
129 | //@: Mean Focal distances [cm]
|
---|
130 | static float ct_Focal_mean;
|
---|
131 |
|
---|
132 | //@: STDev. Focal distances [cm]
|
---|
133 | static float ct_Focal_std;
|
---|
134 |
|
---|
135 | //@: Mean Point Spread function [cm]
|
---|
136 | static float ct_PSpread_mean;
|
---|
137 |
|
---|
138 | //@: STDev. Point Spread function [cm]
|
---|
139 | static float ct_PSpread_std;
|
---|
140 |
|
---|
141 | //@: STDev. Adjustmente deviation [cm]
|
---|
142 | static float ct_Adjustment_std;
|
---|
143 |
|
---|
144 | //@: Radius of the Black Spot in mirror [cm]
|
---|
145 | static float ct_BlackSpot_rad;
|
---|
146 |
|
---|
147 | //@: Radius of one mirror [cm]
|
---|
148 | static float ct_RMirror;
|
---|
149 |
|
---|
150 | //@: Camera width [cm]
|
---|
151 | static float ct_CameraWidth;
|
---|
152 |
|
---|
153 | //@: Pixel width [cm]
|
---|
154 | static float ct_PixelWidth;
|
---|
155 |
|
---|
156 | //@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)
|
---|
157 | static float ct_PixelWidth_corner_2_corner;
|
---|
158 |
|
---|
159 | //@: ct_PixelWidth_corner_2_corner / 2
|
---|
160 | static float ct_PixelWidth_corner_2_corner_half;
|
---|
161 |
|
---|
162 | //@: Number of mirrors
|
---|
163 | static int ct_NMirrors = 0;
|
---|
164 |
|
---|
165 | //@: Number of pixels
|
---|
166 | static int ct_NPixels;
|
---|
167 |
|
---|
168 | //@: Number of pixels
|
---|
169 | static int ct_NCentralPixels;
|
---|
170 |
|
---|
171 | //@: Number of pixels
|
---|
172 | static int ct_NGapPixels;
|
---|
173 |
|
---|
174 | //@: ct_Apot = ct_PixelWidth / 2
|
---|
175 | static float ct_Apot;
|
---|
176 |
|
---|
177 | //@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth
|
---|
178 | static float ct_2Apot;
|
---|
179 |
|
---|
180 | //@: name of the CT definition file to use
|
---|
181 | static char ct_filename[256];
|
---|
182 |
|
---|
183 | //@: list of showers to be skipped
|
---|
184 | static int *Skip;
|
---|
185 |
|
---|
186 | //@: number of showers to be skipped
|
---|
187 | static int nSkip=0;
|
---|
188 |
|
---|
189 | //@: flag: TRUE: data come from STDIN; FALSE: from file
|
---|
190 | static int Data_From_STDIN = FALSE;
|
---|
191 |
|
---|
192 | //@: flag: TRUE: write all images to output; FALSE: only triggered showers
|
---|
193 | static int Write_All_Images = FALSE;
|
---|
194 |
|
---|
195 | //@: flag: TRUE: write all data to output; FALSE: only triggered showers
|
---|
196 | static int Write_All_Data = FALSE;
|
---|
197 |
|
---|
198 | //@: flag: TRUE: selection on the energy
|
---|
199 | static int Select_Energy = TRUE;
|
---|
200 |
|
---|
201 | //@: Lower edge of the selected energy range (in GeV)
|
---|
202 | static float Select_Energy_le = 0.0;
|
---|
203 |
|
---|
204 | //@: Upper edge of the selected energy range (in GeV)
|
---|
205 | static float Select_Energy_ue = 100000.0;
|
---|
206 |
|
---|
207 | //!@}
|
---|
208 |
|
---|
209 | /*!@"
|
---|
210 |
|
---|
211 | The following double-pointer is a 2-dimensional table with information
|
---|
212 | about each pixel. The routine read_pixels will generate
|
---|
213 | the information for filling it using igen_pixel_coordinates().
|
---|
214 |
|
---|
215 | @"*/
|
---|
216 |
|
---|
217 | //!@{
|
---|
218 | // Pointer to a tables/Arrays with information about the pixels
|
---|
219 | // and data stored on them with information about the pixels
|
---|
220 |
|
---|
221 |
|
---|
222 | //@: coordinates x,y for each pixel
|
---|
223 | static float **pixary;
|
---|
224 |
|
---|
225 | //@: indexes of pixels neighbours of a given one
|
---|
226 | static int **pixneig;
|
---|
227 |
|
---|
228 | //@: number of neighbours a pixel have
|
---|
229 | static int *npixneig;
|
---|
230 |
|
---|
231 | //@: contents of the pixels (ph.e.)
|
---|
232 | static float *fnpix;
|
---|
233 |
|
---|
234 | //@: contents of the pixels (ph.e.) after cleanning
|
---|
235 | static float *fnpixclean;
|
---|
236 |
|
---|
237 |
|
---|
238 | //!@}
|
---|
239 |
|
---|
240 | /*!@"
|
---|
241 |
|
---|
242 | The following double-pointer is a 2-dimensional table with the
|
---|
243 | Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
|
---|
244 | of the wavelength @$\lambda@$. The routine |read_pixels()| will read
|
---|
245 | also this information from the file |qe.dat|.
|
---|
246 |
|
---|
247 | @"*/
|
---|
248 |
|
---|
249 | //!@{
|
---|
250 | // Pointer to a table with QE, number of datapoints, and wavelengths
|
---|
251 |
|
---|
252 | //@: table of QE
|
---|
253 | static float ***QE;
|
---|
254 |
|
---|
255 | //@: number of datapoints for the QE curve
|
---|
256 | static int pointsQE;
|
---|
257 |
|
---|
258 | //@: table of QE
|
---|
259 | static float *QElambda;
|
---|
260 | //!@}
|
---|
261 |
|
---|
262 | /*!@"
|
---|
263 |
|
---|
264 | The following double-pointer is a 2-dimensional table with information
|
---|
265 | about each mirror in the dish. The routine |read_ct_file()| will read
|
---|
266 | this information from the CT definition file.
|
---|
267 |
|
---|
268 | @"*/
|
---|
269 |
|
---|
270 | //!@{
|
---|
271 | // Pointer to a table with the following info.:
|
---|
272 |
|
---|
273 | static float **ct_data;
|
---|
274 |
|
---|
275 | /*
|
---|
276 | * TYPE=0 (CT1)
|
---|
277 | * i s rho theta x y z thetan phin xn yn zn
|
---|
278 | *
|
---|
279 | * i : number of the mirror
|
---|
280 | * s : arc length [cm]
|
---|
281 | * rho : polar rho of the position of the center of the mirror [cm]
|
---|
282 | * theta : polar angle of the position of the center of the mirror [cm]
|
---|
283 | * x : x coordinate of the center of the mirror [cm]
|
---|
284 | * y : y coordinate of the center of the mirror [cm]
|
---|
285 | * z : z coordinate of the center of the mirror [cm]
|
---|
286 | * thetan : polar theta angle of the direction where the mirror points to
|
---|
287 | * phin : polar phi angle of the direction where the mirror points to
|
---|
288 | * xn : xn coordinate of the normal vector in the center (normalized)
|
---|
289 | * yn : yn coordinate of the normal vector in the center (normalized)
|
---|
290 | * zn : zn coordinate of the normal vector in the center (normalized)
|
---|
291 | *
|
---|
292 | * TYPE=1 (MAGIC)
|
---|
293 | * i f sx sy x y z thetan phin
|
---|
294 | *
|
---|
295 | * i : number of the mirror
|
---|
296 | * f : focal distance of that mirror
|
---|
297 | * sx : curvilinear coordinate of mirror's center in X[cm]
|
---|
298 | * sy : curvilinear coordinate of mirror's center in X[cm]
|
---|
299 | * x : x coordinate of the center of the mirror [cm]
|
---|
300 | * y : y coordinate of the center of the mirror [cm]
|
---|
301 | * z : z coordinate of the center of the mirror [cm]
|
---|
302 | * thetan : polar theta angle of the direction where the mirror points to
|
---|
303 | * phin : polar phi angle of the direction where the mirror points to
|
---|
304 | * xn : xn coordinate of the normal vector in the center (normalized)
|
---|
305 | * yn : yn coordinate of the normal vector in the center (normalized)
|
---|
306 | * zn : zn coordinate of the normal vector in the center (normalized)
|
---|
307 | */
|
---|
308 | //!@}
|
---|
309 |
|
---|
310 | /*!@"
|
---|
311 |
|
---|
312 | We define a table into where random numbers will be stored.
|
---|
313 | The routines used for random number generation are provided by
|
---|
314 | |RANLIB| (taken from NETLIB, |www.netlib.org|), and by
|
---|
315 | the routine |double drand48(void)| (prototype defined in
|
---|
316 | |stdlib.h|) through the macro |RandomNumber| defined in
|
---|
317 | |camera.h|.
|
---|
318 |
|
---|
319 | @"*/
|
---|
320 |
|
---|
321 | /*!@"
|
---|
322 |
|
---|
323 | The following are the set of parameters calculated for each image.
|
---|
324 | The routines for their calculations are in |moments.cxx|.
|
---|
325 |
|
---|
326 | @"*/
|
---|
327 |
|
---|
328 | //!@{
|
---|
329 | // parameters of the images
|
---|
330 |
|
---|
331 | static Moments_Info *moments_ptr;
|
---|
332 | static LenWid_Info *lenwid_ptr;
|
---|
333 |
|
---|
334 | static float *maxs;
|
---|
335 | static int *nmaxs;
|
---|
336 | static float length, width, dist, xdist, azw, miss, alpha, *conc;
|
---|
337 | static float phiasym, asymx, asymy;
|
---|
338 | static float charge, smax, maxtrigthr_phe;
|
---|
339 |
|
---|
340 | //!@}
|
---|
341 |
|
---|
342 | extern char FileName[];
|
---|
343 |
|
---|
344 |
|
---|
345 | //=-----------------------------------------------------------
|
---|
346 | // @subsection Main program.
|
---|
347 |
|
---|
348 | //!@{
|
---|
349 |
|
---|
350 | //++++++++++++++++++++++++++++++++++++++++
|
---|
351 | // MAIN PROGRAM
|
---|
352 | //----------------------------------------
|
---|
353 |
|
---|
354 | int main(int argc, char **argv)
|
---|
355 | {
|
---|
356 |
|
---|
357 | //!@' @#### Definition of variables.
|
---|
358 | //@'
|
---|
359 |
|
---|
360 | char inname[256]; //@< input file name
|
---|
361 | char starfieldname[256]; //@< starfield input file name
|
---|
362 | char datname[256]; //@< data (ASCII) output file name
|
---|
363 | char diagname[256]; //@< diagnistic output file (ROOT format)
|
---|
364 | char rootname[256] ; //@< ROOT file name
|
---|
365 |
|
---|
366 | char parname[256]; //@< parameters file name
|
---|
367 |
|
---|
368 | char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
|
---|
369 |
|
---|
370 | FILE *inputfile; //@< stream for the input file
|
---|
371 | ofstream datafile; //@< stream for the data file
|
---|
372 |
|
---|
373 | MCEventHeader mcevth; //@< Event Header class (MC)
|
---|
374 |
|
---|
375 | Photoelectron *photoe = NULL; //@< array of the photoelectrons of one event
|
---|
376 | int inumphe; //@< number of photoelectrons in an event
|
---|
377 |
|
---|
378 | float arrtmin_ns; //@ arrival time of the first photoelectron
|
---|
379 | float arrtmax_ns; //@ arrival time of the last photoelectron
|
---|
380 |
|
---|
381 | float thetaCT, phiCT; //@< parameters of a given shower
|
---|
382 | float thetashw, phishw; //@< parameters of a given shower
|
---|
383 | float coreD, coreX, coreY; //@< core position and distance
|
---|
384 | float impactD; //@< impact parameter
|
---|
385 | float l1, m1, n1; //@< auxiliary variables
|
---|
386 | float l2, m2, n2; //@< auxiliary variables
|
---|
387 | float num, den; //@< auxiliary variables
|
---|
388 |
|
---|
389 | int nshow=0; //@< partial number of shower in a given run
|
---|
390 | int ntshow=0; //@< total number of showers
|
---|
391 | int ncph=0; //@< partial number of photons in a given run
|
---|
392 | int ntcph=0; //@< total number of photons
|
---|
393 |
|
---|
394 | int i, j, k; //@< simple counters
|
---|
395 |
|
---|
396 | int simulateNSB; //@< Will we simulate NSB?
|
---|
397 | float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel)
|
---|
398 | float diffnsb_phepns[iMAXNUMPIX]; //@< diffuse NSB values for each pixel derived
|
---|
399 | //@< from meanNSB
|
---|
400 | float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS]; //@< non-diffuse nsb
|
---|
401 | //@< photoelectron rates
|
---|
402 | float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
|
---|
403 | EXTWAVEBAND1,
|
---|
404 | EXTWAVEBAND2,
|
---|
405 | EXTWAVEBAND3,
|
---|
406 | EXTWAVEBAND4,
|
---|
407 | EXTWAVEBAND5
|
---|
408 | };
|
---|
409 | float baseline_mv[iMAXNUMPIX]; //@< The baseline (mV) caused by the NSB; to be subtracted
|
---|
410 | //@< in order to simulate the preamps' AC coupling
|
---|
411 |
|
---|
412 | float qThreshold; //@< Threshold value
|
---|
413 | float qTailCut; //@< Tail Cut value
|
---|
414 | int nIslandsCut; //@< Islands Cut value
|
---|
415 | int countIslands; //@< Will we count the islands?
|
---|
416 | int anaPixels;
|
---|
417 |
|
---|
418 | float fCorrection; //@< Factor to apply to pixel values (def. 1.)
|
---|
419 |
|
---|
420 | int trigger; //@< trigger flag
|
---|
421 | int itrigger; //@< index of pixel fired
|
---|
422 | int ntrigger = 0; //@< number of triggers in the whole file
|
---|
423 | // MTrigger Trigger ; //@< A instance of the Class MTrigger
|
---|
424 |
|
---|
425 | // MMcTrig *McTrig = new MMcTrig() ;
|
---|
426 |
|
---|
427 | float plateScale_cm2deg; //@< plate scale (deg/cm)
|
---|
428 | float degTriggerZone; //@< trigger area in the camera (radius, in deg.)
|
---|
429 |
|
---|
430 | float dtheta, dphi; //@< deviations of CT from shower axis
|
---|
431 |
|
---|
432 | int still_in_loop = FALSE;
|
---|
433 |
|
---|
434 | float *image_data;
|
---|
435 | int nvar, hidt;
|
---|
436 |
|
---|
437 | struct camera cam; // structure holding the camera definition
|
---|
438 |
|
---|
439 |
|
---|
440 | //!@' @#### Definition of variables for |getopt()|.
|
---|
441 | //@'
|
---|
442 |
|
---|
443 | int ch, errflg = 0; //@< used by getopt
|
---|
444 |
|
---|
445 | /*!@'
|
---|
446 |
|
---|
447 | @#### Beginning of the program.
|
---|
448 |
|
---|
449 | We start with the main program. First we (could) make some
|
---|
450 | presentation, and follows the reading of the parameters file (now
|
---|
451 | from the |stdin|), the reading of the CT parameters file, and the
|
---|
452 | creation of the output file, where the processed data will be
|
---|
453 | stored.
|
---|
454 |
|
---|
455 | */
|
---|
456 |
|
---|
457 | //++
|
---|
458 | // START
|
---|
459 | //--
|
---|
460 |
|
---|
461 | // make unbuffered output
|
---|
462 |
|
---|
463 | cout.setf ( ios::stdio );
|
---|
464 |
|
---|
465 | // parse command line options (see reflector.h)
|
---|
466 |
|
---|
467 | parname[0] = '\0';
|
---|
468 |
|
---|
469 | optarg = NULL;
|
---|
470 | while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
|
---|
471 | switch (ch) {
|
---|
472 | case 'f':
|
---|
473 | strcpy(parname, optarg);
|
---|
474 | break;
|
---|
475 | case 'h':
|
---|
476 | usage();
|
---|
477 | break;
|
---|
478 | default :
|
---|
479 | errflg++;
|
---|
480 | }
|
---|
481 |
|
---|
482 | // show help if error
|
---|
483 |
|
---|
484 | if ( errflg>0 )
|
---|
485 | usage();
|
---|
486 |
|
---|
487 | // make some sort of presentation
|
---|
488 |
|
---|
489 | present();
|
---|
490 |
|
---|
491 | // read parameters file
|
---|
492 |
|
---|
493 | if ( strlen(parname) < 1 )
|
---|
494 | readparam(NULL);
|
---|
495 | else
|
---|
496 | readparam(parname);
|
---|
497 |
|
---|
498 | // read data from file or from STDIN?
|
---|
499 |
|
---|
500 | Data_From_STDIN = get_data_from_stdin();
|
---|
501 |
|
---|
502 | // write all images, even those without trigger?
|
---|
503 |
|
---|
504 | Write_All_Images = get_write_all_images();
|
---|
505 |
|
---|
506 | // write all data (i.e., ph.e.s in pixels)
|
---|
507 |
|
---|
508 | Write_All_Data = get_write_all_data();
|
---|
509 |
|
---|
510 | // get filenames
|
---|
511 |
|
---|
512 | strcpy( inname, get_input_filename() );
|
---|
513 | strcpy( starfieldname, get_starfield_filename() );
|
---|
514 | strcpy( datname, get_data_filename() );
|
---|
515 | strcpy( diagname, get_diag_filename() );
|
---|
516 | strcpy( rootname, get_root_filename() );
|
---|
517 | strcpy( ct_filename, get_ct_filename() );
|
---|
518 |
|
---|
519 | // get different parameters of the simulation
|
---|
520 |
|
---|
521 | qThreshold = get_threshold();
|
---|
522 | qTailCut = get_tail_cut();
|
---|
523 | simulateNSB = get_nsb( &meanNSB );
|
---|
524 | countIslands = get_islands_cut( &nIslandsCut );
|
---|
525 |
|
---|
526 | // get selections on the parameters
|
---|
527 |
|
---|
528 | Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
|
---|
529 |
|
---|
530 | // log filenames information
|
---|
531 |
|
---|
532 | log(SIGNATURE,
|
---|
533 | "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n",
|
---|
534 | "Filenames",
|
---|
535 | "In", inname,
|
---|
536 | "Stars", starfieldname,
|
---|
537 | "CT", ct_filename,
|
---|
538 | "Data", datname,
|
---|
539 | "Diag", diagname,
|
---|
540 | "ROOT", rootname
|
---|
541 | );
|
---|
542 |
|
---|
543 |
|
---|
544 | // log flags information
|
---|
545 |
|
---|
546 | log(SIGNATURE,
|
---|
547 | "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
|
---|
548 | "Flags",
|
---|
549 | "Data_From_STDIN", ONoff(Data_From_STDIN),
|
---|
550 | "Write_All_Images", ONoff(Write_All_Images),
|
---|
551 | "Write_All_Data", ONoff(Write_All_Data));
|
---|
552 |
|
---|
553 | // log parameters information
|
---|
554 |
|
---|
555 | log(SIGNATURE,
|
---|
556 | "%s:\n\t%20s: %f\n\t%20s: %f\n\t%20s: %f %s\n\t%20s: %f %s\n",
|
---|
557 | "Parameters",
|
---|
558 | "q0 (Threshold)", qThreshold,
|
---|
559 | "t0 (Tail-cut)", qTailCut,
|
---|
560 | "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB),
|
---|
561 | "i0 (Islands-cut)", nIslandsCut, ONoff(countIslands));
|
---|
562 |
|
---|
563 | // log selections
|
---|
564 |
|
---|
565 | log(SIGNATURE,
|
---|
566 | "%s:\n\t%20s: %s (%f:%f)\n",
|
---|
567 | "Selections:",
|
---|
568 | "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
|
---|
569 |
|
---|
570 | // set all random numbers seeds
|
---|
571 |
|
---|
572 | setall( get_seeds(0), get_seeds(1) );
|
---|
573 |
|
---|
574 | // get list of showers to evt. skip
|
---|
575 |
|
---|
576 | nSkip = get_nskip_showers();
|
---|
577 |
|
---|
578 | if (nSkip > 0) {
|
---|
579 | Skip = new int[ nSkip ];
|
---|
580 | get_skip_showers( Skip );
|
---|
581 |
|
---|
582 | log(SIGNATURE, "There are some showers to skip:\n");
|
---|
583 | for (i=0; i<nSkip; ++i)
|
---|
584 | log(SIGNATURE, "\tshower # %d\n", Skip[i]);
|
---|
585 | }
|
---|
586 |
|
---|
587 | // read parameters from the ct.def file
|
---|
588 |
|
---|
589 | read_ct_file();
|
---|
590 |
|
---|
591 | // read camera setup
|
---|
592 |
|
---|
593 | read_pixels(&cam);
|
---|
594 |
|
---|
595 | // allocate memory for the photoelectrons
|
---|
596 |
|
---|
597 | photoe = new Photoelectron[iMAXNUMPHE];
|
---|
598 |
|
---|
599 | // initialise ROOT
|
---|
600 |
|
---|
601 | TROOT simple("simple", "MAGIC Telescope Monte Carlo");
|
---|
602 |
|
---|
603 | // prepare ROOT tree for the diagnostic data
|
---|
604 |
|
---|
605 | TFile *hfile;
|
---|
606 |
|
---|
607 | hfile = new TFile( diagname,"RECREATE", "MAGIC Telescope MC diagnostic data");
|
---|
608 |
|
---|
609 | // Create the ROOT Tree for the diagnostic data
|
---|
610 |
|
---|
611 | TTree *tree = new TTree("T","MAGIC Telescope MC diagnostic data");
|
---|
612 | tree->SetAutoSave(100000000);
|
---|
613 |
|
---|
614 | Int_t split = 1;
|
---|
615 | Int_t bsize = 64000;
|
---|
616 | MDiagEventobject *event = 0;
|
---|
617 |
|
---|
618 | // Create one branch. If splitlevel is set, event is a superbranch
|
---|
619 | // creating a sub branch for each data member of the Eventobject event.
|
---|
620 |
|
---|
621 | tree->Branch("event", "MDiagEventobject", &event, bsize, split);
|
---|
622 |
|
---|
623 | // Prepare the raw data output
|
---|
624 |
|
---|
625 | MRawEvt *Evt = new MRawEvt() ;
|
---|
626 | MMcEvt *McEvt = new MMcEvt ();
|
---|
627 |
|
---|
628 | // initalize the ROOT file
|
---|
629 | //
|
---|
630 |
|
---|
631 | TFile outfile ( rootname , "RECREATE" );
|
---|
632 |
|
---|
633 | // create a Tree for the Event data stream
|
---|
634 |
|
---|
635 | TTree EvtTree("EvtTree","Events of Run");
|
---|
636 |
|
---|
637 | bsize=128000; split=1;
|
---|
638 |
|
---|
639 | EvtTree.Branch("MRawEvt","MRawEvt",
|
---|
640 | &Evt, bsize, split);
|
---|
641 |
|
---|
642 | EvtTree.Branch("MMcEvt","MMcEvt",
|
---|
643 | &McEvt, bsize, split);
|
---|
644 |
|
---|
645 |
|
---|
646 | // for safety and for dimensioning image_data: count the elements in the
|
---|
647 | // diagnostic data branch
|
---|
648 |
|
---|
649 | i=0;
|
---|
650 | i++; // "n"
|
---|
651 | i++; // "primary"
|
---|
652 | i++; // "energy"
|
---|
653 | i++; // "cored"
|
---|
654 | i++; // "impact"
|
---|
655 | i++; // "xcore"
|
---|
656 | i++; // "ycore"
|
---|
657 | i++; // "theta"
|
---|
658 | i++; // "phi"
|
---|
659 | i++; // "deviations"
|
---|
660 | i++; // "dtheta"
|
---|
661 | i++; // "dphi"
|
---|
662 | i++; // "trigger"
|
---|
663 | i++; // "ncphs"
|
---|
664 | i++; // "maxpassthr_phe"
|
---|
665 | i++; // "nphes"
|
---|
666 | i++; // "nphes2"
|
---|
667 | i++; // "length"
|
---|
668 | i++; // "width"
|
---|
669 | i++; // "dist"
|
---|
670 | i++; // "xdist"
|
---|
671 | i++; // "azw"
|
---|
672 | i++; // "miss"
|
---|
673 | i++; // "alpha"
|
---|
674 | i++; // "conc2"
|
---|
675 | i++; // "conc3"
|
---|
676 | i++; // "conc4"
|
---|
677 | i++; // "conc5"
|
---|
678 | i++; // "conc6"
|
---|
679 | i++; // "conc7"
|
---|
680 | i++; // "conc8"
|
---|
681 | i++; // "conc9"
|
---|
682 | i++; // "conc10"
|
---|
683 | i++; // "asymx"
|
---|
684 | i++; // "asymy"
|
---|
685 | i++; // "phiasym"
|
---|
686 |
|
---|
687 | nvar = i;
|
---|
688 | image_data = new float[nvar];
|
---|
689 |
|
---|
690 | // set plate scale (deg/cm) and trigger area (deg)
|
---|
691 |
|
---|
692 | plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;
|
---|
693 |
|
---|
694 | if ( ! get_trigger_radius( °TriggerZone ) )
|
---|
695 | degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
|
---|
696 |
|
---|
697 | if ( ! get_correction( &fCorrection ) )
|
---|
698 | fCorrection = 1.0;
|
---|
699 |
|
---|
700 | // number of pixels for parameters
|
---|
701 |
|
---|
702 | anaPixels = get_ana_pixels();
|
---|
703 | anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
|
---|
704 |
|
---|
705 | // prepare the NSB simulation
|
---|
706 |
|
---|
707 | if( simulateNSB ){
|
---|
708 |
|
---|
709 | //
|
---|
710 | // Calculate the non-diffuse NSB photoelectron rates
|
---|
711 | //
|
---|
712 |
|
---|
713 | k = produce_nsbrates( starfieldname,
|
---|
714 | &cam,
|
---|
715 | photoe, // only a dummy here
|
---|
716 | nsbrate_phepns );
|
---|
717 | if (k != 0){
|
---|
718 | cout << "Error when reading starfield... \nExiting.\n";
|
---|
719 | exit(1);
|
---|
720 | }
|
---|
721 |
|
---|
722 | // for(i=0; i<cam.inumpixels; i++){
|
---|
723 | // cout << i;
|
---|
724 | // for(j=0; j<iNUMWAVEBANDS; j++){
|
---|
725 | // cout << " " << j << " " << nsbrate_phepns[i][j];
|
---|
726 | // }
|
---|
727 | // cout << "\n";
|
---|
728 | // }
|
---|
729 |
|
---|
730 | // calculate diffuse rate correcting for the pixel size
|
---|
731 |
|
---|
732 | for(i=0; i<cam.inumpixels; i++){
|
---|
733 | diffnsb_phepns[i] = meanNSB *
|
---|
734 | cam.dpixsizefactor[i] * cam.dpixsizefactor[i];
|
---|
735 | }
|
---|
736 |
|
---|
737 | }
|
---|
738 |
|
---|
739 | //
|
---|
740 | // Read the reflector file with the Cherenkov data
|
---|
741 | //
|
---|
742 |
|
---|
743 | // select input file
|
---|
744 |
|
---|
745 | if ( Data_From_STDIN ) {
|
---|
746 |
|
---|
747 | inputfile = stdin;
|
---|
748 |
|
---|
749 | }
|
---|
750 | else{
|
---|
751 |
|
---|
752 | log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
|
---|
753 | inputfile = fopen( inname, "r" );
|
---|
754 | if ( inputfile == NULL )
|
---|
755 | error( SIGNATURE, "Cannot open input file: %s\n", inname );
|
---|
756 |
|
---|
757 | }
|
---|
758 |
|
---|
759 | // get signature, and check it
|
---|
760 |
|
---|
761 | if(check_reflector_file( inputfile )==FALSE){
|
---|
762 | exit(1);
|
---|
763 | }
|
---|
764 |
|
---|
765 | // open data file
|
---|
766 |
|
---|
767 | log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
|
---|
768 | datafile.open( datname );
|
---|
769 |
|
---|
770 | if ( datafile.bad() )
|
---|
771 | error( SIGNATURE, "Cannot open data file: %s\n", datname );
|
---|
772 |
|
---|
773 | // initializes flag
|
---|
774 |
|
---|
775 | strcpy( flag, " \0" );
|
---|
776 |
|
---|
777 | // allocate space for PMTs numbers of pixels
|
---|
778 |
|
---|
779 | fnpix = new float [ ct_NPixels ];
|
---|
780 | fnpixclean = new float [ ct_NPixels ];
|
---|
781 |
|
---|
782 | moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
|
---|
783 |
|
---|
784 | // initialize baseline
|
---|
785 |
|
---|
786 | for(i=0; i<cam.inumpixels; i++){
|
---|
787 | baseline_mv[i] = 0.;
|
---|
788 | }
|
---|
789 |
|
---|
790 | // Instance of the Trigger class
|
---|
791 |
|
---|
792 | MTrigger Trigger ; //@< A instance of the Class MTrigger
|
---|
793 |
|
---|
794 | //!@' @#### Main loop.
|
---|
795 | //@'
|
---|
796 |
|
---|
797 | // get flag
|
---|
798 |
|
---|
799 | fread( flag, SIZE_OF_FLAGS, 1, inputfile );
|
---|
800 |
|
---|
801 | // loop over the file
|
---|
802 |
|
---|
803 | still_in_loop = TRUE;
|
---|
804 |
|
---|
805 | while (
|
---|
806 | ((! Data_From_STDIN) && ( !feof(inputfile) ))
|
---|
807 | ||
|
---|
808 | (Data_From_STDIN && still_in_loop)
|
---|
809 | ) {
|
---|
810 |
|
---|
811 | // reading .rfl files
|
---|
812 | if(!isA( flag, FLAG_START_OF_RUN )){
|
---|
813 | error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
|
---|
814 | }
|
---|
815 | else { // found start of run
|
---|
816 |
|
---|
817 | nshow=0;
|
---|
818 |
|
---|
819 | fread( flag, SIZE_OF_FLAGS, 1, inputfile );
|
---|
820 |
|
---|
821 | while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
|
---|
822 |
|
---|
823 | ++nshow;
|
---|
824 | log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
|
---|
825 |
|
---|
826 | // get MCEventHeader
|
---|
827 |
|
---|
828 | fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
|
---|
829 |
|
---|
830 | // calculate core distance and impact parameter
|
---|
831 |
|
---|
832 | coreD = mcevth.get_core(&coreX, &coreY);
|
---|
833 |
|
---|
834 | // calculate impact parameter (shortest distance betwee the original
|
---|
835 | // trajectory of the primary (assumed shower-axis) and the
|
---|
836 | // direction where the telescope points to
|
---|
837 | //
|
---|
838 | // we use the following equation, given that the shower core position
|
---|
839 | // is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
|
---|
840 | // and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
|
---|
841 | // and (l2,m2,n2)
|
---|
842 | //
|
---|
843 | // | |
|
---|
844 | // | x1-x2 y1-y2 z1-z2 |
|
---|
845 | // | |
|
---|
846 | // + | l1 m1 n1 |
|
---|
847 | // - | |
|
---|
848 | // | l2 m2 n2 |
|
---|
849 | // | |
|
---|
850 | // dist = ------------------------------------ ( > 0 )
|
---|
851 | // [ |l1 m1|2 |m1 n1|2 |n1 l1|2 ]1/2
|
---|
852 | // [ | | + | | + | | ]
|
---|
853 | // [ |l2 m2| |m2 n2| |n2 l2| ]
|
---|
854 | //
|
---|
855 | // playing a little bit, we get this reduced for in our case:
|
---|
856 | //
|
---|
857 | //
|
---|
858 | // dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
|
---|
859 | // [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
|
---|
860 | // 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
|
---|
861 |
|
---|
862 | // read the direction of the incoming shower
|
---|
863 |
|
---|
864 | thetashw = mcevth.get_theta();
|
---|
865 | phishw = mcevth.get_phi();
|
---|
866 |
|
---|
867 | // calculate vector for shower
|
---|
868 |
|
---|
869 | l1 = sin(thetashw)*cos(phishw);
|
---|
870 | m1 = sin(thetashw)*sin(phishw);
|
---|
871 | n1 = cos(thetashw);
|
---|
872 |
|
---|
873 | // read the deviation of the telescope with respect to the shower
|
---|
874 |
|
---|
875 | mcevth.get_deviations ( &thetaCT, &phiCT );
|
---|
876 |
|
---|
877 | if ( (thetaCT == 0.) && (phiCT == 0.) ) {
|
---|
878 |
|
---|
879 | // CT was looking to the source (both lines are parallel)
|
---|
880 | // therefore, we calculate the impact parameter as the distance
|
---|
881 | // between the CT axis and the core position
|
---|
882 |
|
---|
883 | impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
|
---|
884 |
|
---|
885 | } else {
|
---|
886 |
|
---|
887 | // the shower comes off-axis
|
---|
888 |
|
---|
889 | // obtain with this the final direction of the CT
|
---|
890 |
|
---|
891 | thetaCT += thetashw;
|
---|
892 | phiCT += phishw;
|
---|
893 |
|
---|
894 | // calculate vector for telescope
|
---|
895 |
|
---|
896 | l2 = sin(thetaCT)*cos(phiCT);
|
---|
897 | m2 = sin(thetaCT)*sin(phiCT);
|
---|
898 | n2 = cos(thetaCT);
|
---|
899 |
|
---|
900 | num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
|
---|
901 | den = (SQR(l1*m2 - l2*m1) +
|
---|
902 | SQR(m1*n2 - m2*n1) +
|
---|
903 | SQR(n1*l2 - n2*l1));
|
---|
904 | den = sqrt(den);
|
---|
905 |
|
---|
906 | impactD = fabs(num)/den;
|
---|
907 |
|
---|
908 | // fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ",
|
---|
909 | // thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2,
|
---|
910 | // num, den);
|
---|
911 |
|
---|
912 | }
|
---|
913 |
|
---|
914 | // read the photons and produce the photoelectrons
|
---|
915 |
|
---|
916 | k = produce_phes( inputfile,
|
---|
917 | &cam,
|
---|
918 | WAVEBANDBOUND1,
|
---|
919 | WAVEBANDBOUND6,
|
---|
920 | photoe, // will be changed by the function!
|
---|
921 | &inumphe, // important for later: the size of photoe[]
|
---|
922 | fnpix, // will be changed by the function!
|
---|
923 | &ncph, // will be changed by the function!
|
---|
924 | &arrtmin_ns, // will be changed by the function!
|
---|
925 | &arrtmax_ns // will be changed by the function!
|
---|
926 | );
|
---|
927 |
|
---|
928 | if( k != 0 ){ // non-zero returnvalue means error
|
---|
929 | cout << "Exiting.\n";
|
---|
930 | exit(1);
|
---|
931 | }
|
---|
932 |
|
---|
933 | log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
|
---|
934 | ncph, ntcph);
|
---|
935 |
|
---|
936 | ntcph += ncph;
|
---|
937 |
|
---|
938 | // skip it ?
|
---|
939 |
|
---|
940 | for ( i=0; i<nSkip; ++i ) {
|
---|
941 | if (Skip[i] == (nshow+ntshow)) {
|
---|
942 | i = -1;
|
---|
943 | break;
|
---|
944 | }
|
---|
945 | }
|
---|
946 |
|
---|
947 | // if after the previous loop, the exit value of i is -1
|
---|
948 | // then the shower number is in the list of showers to be
|
---|
949 | // skipped
|
---|
950 |
|
---|
951 | if (i == -1) {
|
---|
952 | log(SIGNATURE, "\t\tskipped!\n");
|
---|
953 | continue;
|
---|
954 | }
|
---|
955 |
|
---|
956 | // energy cut
|
---|
957 |
|
---|
958 | if ( Select_Energy ) {
|
---|
959 | if (( mcevth.get_energy() < Select_Energy_le ) ||
|
---|
960 | ( mcevth.get_energy() > Select_Energy_ue )) {
|
---|
961 | log(SIGNATURE, "select_energy: shower rejected.\n");
|
---|
962 | continue;
|
---|
963 | }
|
---|
964 | }
|
---|
965 |
|
---|
966 | // NSB simulation
|
---|
967 |
|
---|
968 | if(simulateNSB){
|
---|
969 |
|
---|
970 | k = produce_nsb_phes( &arrtmin_ns, // will be changed by the function!
|
---|
971 | &arrtmax_ns, // will be changed by the function!
|
---|
972 | thetaCT,
|
---|
973 | &cam,
|
---|
974 | nsbrate_phepns,
|
---|
975 | diffnsb_phepns,
|
---|
976 | ext,
|
---|
977 | fnpix, // will be changed by the function!
|
---|
978 | photoe, // will be changed by the function!
|
---|
979 | &inumphe, // important for later: the size of photoe[]
|
---|
980 | baseline_mv // will be generated by the function
|
---|
981 | );
|
---|
982 |
|
---|
983 | if( k != 0 ){ // non-zero returnvalue means error
|
---|
984 | cout << "Exiting.\n";
|
---|
985 | exit(1);
|
---|
986 | }
|
---|
987 |
|
---|
988 | }// end if(simulateNSB) ...
|
---|
989 |
|
---|
990 |
|
---|
991 | // cout << arrtmin_ns << " " << arrtmax_ns << "\n";
|
---|
992 | // for(i=0; i<cam.inumpixels; i++){
|
---|
993 | // cout << i << " " << baseline_mv[i] <<"\n";
|
---|
994 | // }
|
---|
995 |
|
---|
996 | cout << "Total number of phes: " << inumphe <<endl;
|
---|
997 |
|
---|
998 | // TRIGGER HERE
|
---|
999 |
|
---|
1000 | //
|
---|
1001 | // remove the prvious values of the analog signal for
|
---|
1002 | // each pixel and put the new ones from photoe array
|
---|
1003 | //
|
---|
1004 | Trigger.Reset();
|
---|
1005 |
|
---|
1006 | for(i=0;i<inumphe;i++){
|
---|
1007 | Trigger.FillShow(photoe[i].ipixnum,float((photoe[i].iarrtime_ns-arrtmin_ns+TOTAL_TRIGGER_TIME/10.0)));
|
---|
1008 | }
|
---|
1009 |
|
---|
1010 | //
|
---|
1011 | // now the noise of the electronic
|
---|
1012 | // (preamps, optical transmission,..) is introduced.
|
---|
1013 | // This is done inside the class MTrigger by the method ElecNoise.
|
---|
1014 | //
|
---|
1015 | Trigger.ElecNoise() ;
|
---|
1016 |
|
---|
1017 | //
|
---|
1018 | // look if in all the signals in the trigger signal branch
|
---|
1019 | // is a possible Trigger. Therefore we habe to diskriminate all
|
---|
1020 | // the simulated analog signals (Method Diskriminate in class
|
---|
1021 | // MTrigger). We look simultanously for the moments at which
|
---|
1022 | // there are more than TRIGGER_MULTI pixels above the
|
---|
1023 | // CHANNEL_THRESHOLD.
|
---|
1024 | //
|
---|
1025 |
|
---|
1026 | Trigger.Diskriminate();
|
---|
1027 |
|
---|
1028 | //McTrig->SetZeroLevel(trigger = (Short_t) Trigger.ZeroLevel()) ;
|
---|
1029 | trigger = (Short_t) Trigger.ZeroLevel();
|
---|
1030 |
|
---|
1031 | //
|
---|
1032 | // Start the First Level Trigger simulation
|
---|
1033 | //
|
---|
1034 | if ( trigger > 0 ) {
|
---|
1035 | // McTrig->SetFirstLevel ( trigger = Trigger.FirstLevel() ) ;
|
---|
1036 | trigger = Trigger.FirstLevel();
|
---|
1037 | }
|
---|
1038 | if ( trigger>0 ) {
|
---|
1039 | cout << "TRIG"<<" "<<thetashw <<" "<< trigger << endl;
|
---|
1040 | }
|
---|
1041 |
|
---|
1042 | // for( i=0; i<inumphe; i++){
|
---|
1043 | // cout << "phe " << photoe[i].ipixnum << " " << photoe[i].iarrtime_ns << "\n";
|
---|
1044 | // }
|
---|
1045 |
|
---|
1046 | if ( trigger>0 ) {
|
---|
1047 |
|
---|
1048 | itrigger = i;
|
---|
1049 | ++ntrigger;
|
---|
1050 |
|
---|
1051 | memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels );
|
---|
1052 |
|
---|
1053 |
|
---|
1054 | #ifdef __ISLANDS__
|
---|
1055 |
|
---|
1056 | //!@' @#### Islands algorithm.
|
---|
1057 | //@'
|
---|
1058 |
|
---|
1059 | //++
|
---|
1060 | // islands counting, and cleanning
|
---|
1061 | //--
|
---|
1062 |
|
---|
1063 | if ( countIslands )
|
---|
1064 | do_islands( ct_NPixels, fnpixclean, pixneig, npixneig,
|
---|
1065 | countIslands, nIslandsCut);
|
---|
1066 |
|
---|
1067 | #endif // __ISLANDS__
|
---|
1068 |
|
---|
1069 |
|
---|
1070 | //!@' @#### Calculation of parameters of the image.
|
---|
1071 | //@'
|
---|
1072 |
|
---|
1073 | //++
|
---|
1074 | // moments calculation
|
---|
1075 | //--
|
---|
1076 |
|
---|
1077 | // calculate moments and other things
|
---|
1078 |
|
---|
1079 | moments_ptr = moments( anaPixels, fnpixclean, pixary,
|
---|
1080 | plateScale_cm2deg, 0 );
|
---|
1081 |
|
---|
1082 | charge = moments_ptr->charge ;
|
---|
1083 | smax = moments_ptr->smax ;
|
---|
1084 | maxs = moments_ptr->maxs ;
|
---|
1085 | nmaxs = moments_ptr->nmaxs ;
|
---|
1086 | length = moments_ptr->length ;
|
---|
1087 | width = moments_ptr->width ;
|
---|
1088 | dist = moments_ptr->dist ;
|
---|
1089 | xdist = moments_ptr->xdist ;
|
---|
1090 | azw = moments_ptr->azw ;
|
---|
1091 | miss = moments_ptr->miss ;
|
---|
1092 | alpha = moments_ptr->alpha ;
|
---|
1093 | conc = moments_ptr->conc ;
|
---|
1094 | asymx = moments_ptr->asymx ;
|
---|
1095 | asymx = moments_ptr->asymx ;
|
---|
1096 | phiasym= moments_ptr->phi;
|
---|
1097 |
|
---|
1098 | lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary,
|
---|
1099 | plateScale_cm2deg,
|
---|
1100 | ct_PixelWidth_corner_2_corner_half);
|
---|
1101 |
|
---|
1102 |
|
---|
1103 | // fill the diagnostic Tree
|
---|
1104 |
|
---|
1105 | event = new MDiagEventobject();
|
---|
1106 |
|
---|
1107 | i=0;
|
---|
1108 | image_data[i] = event->n = hidt/10; i++;
|
---|
1109 | image_data[i] = event->primary = mcevth.get_primary(); i++;
|
---|
1110 | image_data[i] = event->energy = mcevth.get_energy(); i++;
|
---|
1111 | image_data[i] = event->cored = coreD; i++;
|
---|
1112 | image_data[i] = event->impact = impactD; i++;
|
---|
1113 | image_data[i] = event->xcore = coreX; i++;
|
---|
1114 | image_data[i] = event->ycore = coreY; i++;
|
---|
1115 | image_data[i] = event->theta = mcevth.get_theta(); i++;
|
---|
1116 | image_data[i] = event->phi = mcevth.get_phi(); i++;
|
---|
1117 | image_data[i] = event->deviations = mcevth.get_deviations (&dtheta, &dphi); i++;
|
---|
1118 | image_data[i] = event->dtheta = dtheta; i++;
|
---|
1119 | image_data[i] = event->dphi = dphi; i++;
|
---|
1120 | image_data[i] = event->trigger = trigger; i++;
|
---|
1121 | image_data[i] = event->ncphs = ncph; i++;
|
---|
1122 | image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++;
|
---|
1123 | image_data[i] = event->nphes = charge; i++;
|
---|
1124 | image_data[i] = event->nphes2 = smax; i++;
|
---|
1125 | image_data[i] = event->length = length; i++;
|
---|
1126 | image_data[i] = event->width = width; i++;
|
---|
1127 | image_data[i] = event->dist = dist; i++;
|
---|
1128 | image_data[i] = event->xdist = xdist; i++;
|
---|
1129 | image_data[i] = event->azw = azw; i++;
|
---|
1130 | image_data[i] = event->miss = miss; i++;
|
---|
1131 | image_data[i] = event->alpha = alpha; i++;
|
---|
1132 | image_data[i] = event->conc2 = conc[0]; i++;
|
---|
1133 | image_data[i] = event->conc3 = conc[1]; i++;
|
---|
1134 | image_data[i] = event->conc4 = conc[2]; i++;
|
---|
1135 | image_data[i] = event->conc5 = conc[3]; i++;
|
---|
1136 | image_data[i] = event->conc6 = conc[4]; i++;
|
---|
1137 | image_data[i] = event->conc7 = conc[5]; i++;
|
---|
1138 | image_data[i] = event->conc8 = conc[6]; i++;
|
---|
1139 | image_data[i] = event->conc9 = conc[7]; i++;
|
---|
1140 | image_data[i] = event->conc10 = conc[8]; i++;
|
---|
1141 | image_data[i] = event->asymx = asymx; i++;
|
---|
1142 | image_data[i] = event->asymy = asymy; i++;
|
---|
1143 | image_data[i] = event->phiasym = phiasym; i++;
|
---|
1144 |
|
---|
1145 | // there should be "nvar" variables
|
---|
1146 |
|
---|
1147 | if ( i != nvar )
|
---|
1148 | error( SIGNATURE, "Wrong entry number for diagnostic data.\n" );
|
---|
1149 |
|
---|
1150 | tree->Fill();
|
---|
1151 | delete event;
|
---|
1152 |
|
---|
1153 | // put information in the data file,
|
---|
1154 |
|
---|
1155 | datafile << ntrigger;
|
---|
1156 | for(i=0;i<nvar;i++) {
|
---|
1157 | datafile << ' ' << image_data[i];
|
---|
1158 | }
|
---|
1159 |
|
---|
1160 | // revert the fnpixclean matrix into fnpix
|
---|
1161 | // (now we do this, but maybe in a future we want to
|
---|
1162 | // use both fnpix and fnpixclean for different things
|
---|
1163 |
|
---|
1164 | memcpy( fnpix, fnpixclean, sizeof(float) * ct_NPixels );
|
---|
1165 |
|
---|
1166 | // put this information in the data file,
|
---|
1167 |
|
---|
1168 | if ( Write_All_Data ) {
|
---|
1169 | datafile << ' ' << -9999;
|
---|
1170 | for ( i=0; i<ct_NPixels; ++i )
|
---|
1171 | datafile << ' ' << fnpix[i];
|
---|
1172 | }
|
---|
1173 |
|
---|
1174 | datafile << endl;
|
---|
1175 |
|
---|
1176 | mcevth.set_trigger( TRUE );
|
---|
1177 |
|
---|
1178 | log(SIGNATURE, "TRIGGER\n");
|
---|
1179 |
|
---|
1180 | } else { // ( trigger == FALSE )
|
---|
1181 |
|
---|
1182 | event = new MDiagEventobject();
|
---|
1183 |
|
---|
1184 | i=0;
|
---|
1185 | image_data[i] = event->n = hidt/10; i++;
|
---|
1186 | image_data[i] = event->primary = mcevth.get_primary(); i++;
|
---|
1187 | image_data[i] = event->energy = mcevth.get_energy(); i++;
|
---|
1188 | image_data[i] = event->cored = coreD = mcevth.get_core(&coreX, &coreY); i++;
|
---|
1189 | image_data[i] = event->impact = coreD; i++;
|
---|
1190 | image_data[i] = event->xcore = coreX; i++;
|
---|
1191 | image_data[i] = event->ycore = coreY; i++;
|
---|
1192 | image_data[i] = event->theta = mcevth.get_theta(); i++;
|
---|
1193 | image_data[i] = event->phi = mcevth.get_phi(); i++;
|
---|
1194 | image_data[i] = event->deviations = mcevth.get_deviations(&dtheta, &dphi); i++;
|
---|
1195 | image_data[i] = event->dtheta = dtheta; i++;
|
---|
1196 | image_data[i] = event->dphi = dphi; i++;
|
---|
1197 | image_data[i] = event->trigger = trigger; i++;
|
---|
1198 | image_data[i] = event->ncphs = ncph; i++;
|
---|
1199 | image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++;
|
---|
1200 | image_data[i] = -1.; i++;
|
---|
1201 | image_data[i] = -1.; i++;
|
---|
1202 | image_data[i] = -1.; i++;
|
---|
1203 | image_data[i] = -1.; i++;
|
---|
1204 | image_data[i] = -1.; i++;
|
---|
1205 | image_data[i] = -1.; i++;
|
---|
1206 | image_data[i] = -1.; i++;
|
---|
1207 | image_data[i] = -1.; i++;
|
---|
1208 | image_data[i] = -1.; i++;
|
---|
1209 | image_data[i] = -1.; i++;
|
---|
1210 | image_data[i] = -1.; i++;
|
---|
1211 | image_data[i] = -1.; i++;
|
---|
1212 | image_data[i] = -1.; i++;
|
---|
1213 | image_data[i] = -1.; i++;
|
---|
1214 | image_data[i] = -1.; i++;
|
---|
1215 | image_data[i] = -1.; i++;
|
---|
1216 | image_data[i] = -1.; i++;
|
---|
1217 | image_data[i] = -1.; i++;
|
---|
1218 | image_data[i] = -1.; i++;
|
---|
1219 | image_data[i] = -1.; i++;
|
---|
1220 | image_data[i] = -1.; i++;
|
---|
1221 |
|
---|
1222 | // there should be "nvar" variables
|
---|
1223 |
|
---|
1224 | if ( i != nvar )
|
---|
1225 | error( SIGNATURE, "Wrong entry length for Ntuple.\n" );
|
---|
1226 |
|
---|
1227 | tree->Fill();
|
---|
1228 | delete event;
|
---|
1229 |
|
---|
1230 | // put this information in the data file,
|
---|
1231 |
|
---|
1232 | if ( Write_All_Data ) {
|
---|
1233 |
|
---|
1234 | datafile << ntrigger;
|
---|
1235 | for ( i=0; i<nvar; ++i )
|
---|
1236 | datafile << ' ' << image_data[i];
|
---|
1237 |
|
---|
1238 | datafile << -9999;
|
---|
1239 | for ( i=0; i<ct_NPixels; ++i )
|
---|
1240 | datafile << ' ' << fnpix[i];
|
---|
1241 |
|
---|
1242 | datafile << endl;
|
---|
1243 | }
|
---|
1244 |
|
---|
1245 | mcevth.set_trigger( FALSE );
|
---|
1246 |
|
---|
1247 | } // trigger == FALSE
|
---|
1248 |
|
---|
1249 | //!@' @#### Save data.
|
---|
1250 | //@'
|
---|
1251 |
|
---|
1252 | //++++++++++++++++++++++++++++++++++++++++++++++++++
|
---|
1253 | // we now have all information we want
|
---|
1254 | // the only thing we must do now is writing it to
|
---|
1255 | // the output file
|
---|
1256 | //--------------------------------------------------
|
---|
1257 |
|
---|
1258 | //++
|
---|
1259 | // save the image to the file
|
---|
1260 | //--
|
---|
1261 |
|
---|
1262 |
|
---|
1263 |
|
---|
1264 | // look for the next event
|
---|
1265 |
|
---|
1266 | fread( flag, SIZE_OF_FLAGS, 1, inputfile );
|
---|
1267 |
|
---|
1268 | } // end while there is a next event
|
---|
1269 |
|
---|
1270 | if( !isA( flag, FLAG_END_OF_RUN )){
|
---|
1271 | error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
|
---|
1272 | }
|
---|
1273 | else { // found end of run
|
---|
1274 | ntshow += nshow;
|
---|
1275 | log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
|
---|
1276 |
|
---|
1277 | fread( flag, SIZE_OF_FLAGS, 1, inputfile );
|
---|
1278 |
|
---|
1279 | if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
|
---|
1280 | log(SIGNATURE, "End of file . . .\n");
|
---|
1281 | still_in_loop = FALSE;
|
---|
1282 |
|
---|
1283 | if ((! Data_From_STDIN) && ( !feof(inputfile) )){
|
---|
1284 |
|
---|
1285 | // we have concatenated input files.
|
---|
1286 | // get signature of the next part and check it.
|
---|
1287 |
|
---|
1288 | if(check_reflector_file( inputfile )==FALSE){
|
---|
1289 | exit(1);
|
---|
1290 | }
|
---|
1291 |
|
---|
1292 | }
|
---|
1293 |
|
---|
1294 | } // end if found end of file
|
---|
1295 | } // end if found end of run
|
---|
1296 |
|
---|
1297 | fread( flag, SIZE_OF_FLAGS, 1, inputfile );
|
---|
1298 |
|
---|
1299 | } // end if else found start of run
|
---|
1300 | } // end big while loop
|
---|
1301 |
|
---|
1302 | //++
|
---|
1303 | // put the Event to the root file
|
---|
1304 | //--
|
---|
1305 |
|
---|
1306 | EvtTree.Write() ;
|
---|
1307 | outfile.Write() ;
|
---|
1308 | outfile.Close() ;
|
---|
1309 |
|
---|
1310 |
|
---|
1311 | // close input file
|
---|
1312 |
|
---|
1313 | log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
|
---|
1314 | ntshow, ntcph );
|
---|
1315 | log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
|
---|
1316 | ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
|
---|
1317 |
|
---|
1318 | // close files
|
---|
1319 |
|
---|
1320 | log( SIGNATURE, "Closing files\n" );
|
---|
1321 |
|
---|
1322 | if( ! Data_From_STDIN ){
|
---|
1323 | fclose( inputfile );
|
---|
1324 | }
|
---|
1325 | datafile.close();
|
---|
1326 |
|
---|
1327 | hfile->Write();
|
---|
1328 |
|
---|
1329 | hfile->Close();
|
---|
1330 |
|
---|
1331 | // program finished
|
---|
1332 |
|
---|
1333 | log( SIGNATURE, "Done.\n");
|
---|
1334 |
|
---|
1335 | return( 0 );
|
---|
1336 | }
|
---|
1337 | //!@}
|
---|
1338 |
|
---|
1339 | // @T \newpage
|
---|
1340 |
|
---|
1341 | //!@subsection Functions definition.
|
---|
1342 |
|
---|
1343 | //!-----------------------------------------------------------
|
---|
1344 | // @name present
|
---|
1345 | //
|
---|
1346 | // @desc Make some presentation
|
---|
1347 | //
|
---|
1348 | // @date Sat Jun 27 05:58:56 MET DST 1998
|
---|
1349 | //------------------------------------------------------------
|
---|
1350 | // @function
|
---|
1351 |
|
---|
1352 | //!@{
|
---|
1353 | void
|
---|
1354 | present(void)
|
---|
1355 | {
|
---|
1356 | cout << "##################################################\n"
|
---|
1357 | << SIGNATURE << '\n' << '\n'
|
---|
1358 | << "Processor of the reflector output\n"
|
---|
1359 | << "J C Gonzalez, Jun 1998\n"
|
---|
1360 | << "##################################################\n\n"
|
---|
1361 | << flush ;
|
---|
1362 | }
|
---|
1363 | //!@}
|
---|
1364 |
|
---|
1365 |
|
---|
1366 | //!-----------------------------------------------------------
|
---|
1367 | // @name usage
|
---|
1368 | //
|
---|
1369 | // @desc show help
|
---|
1370 | //
|
---|
1371 | // @date Tue Dec 15 16:23:30 MET 1998
|
---|
1372 | //------------------------------------------------------------
|
---|
1373 | // @function
|
---|
1374 |
|
---|
1375 | //!@{
|
---|
1376 | void
|
---|
1377 | usage(void)
|
---|
1378 | {
|
---|
1379 | present();
|
---|
1380 | cout << "\nusage ::\n\n"
|
---|
1381 | << "\t camera "
|
---|
1382 | << " [ -@ paramfile ] "
|
---|
1383 | << " [ -h ] "
|
---|
1384 | << "\n\n or \n\n"
|
---|
1385 | << "\t camera < paramfile"
|
---|
1386 | << "\n\n";
|
---|
1387 | exit(0);
|
---|
1388 | }
|
---|
1389 | //!@}
|
---|
1390 |
|
---|
1391 |
|
---|
1392 | //!-----------------------------------------------------------
|
---|
1393 | // @name log
|
---|
1394 | //
|
---|
1395 | // @desc function to send log information
|
---|
1396 | //
|
---|
1397 | // @var funct Name of the caller function
|
---|
1398 | // @var fmt Format to be used (message)
|
---|
1399 | // @var ... Other information to be shown
|
---|
1400 | //
|
---|
1401 | // @date Sat Jun 27 05:58:56 MET DST 1998
|
---|
1402 | //------------------------------------------------------------
|
---|
1403 | // @function
|
---|
1404 |
|
---|
1405 | //!@{
|
---|
1406 | void
|
---|
1407 | log(const char *funct, char *fmt, ...)
|
---|
1408 | {
|
---|
1409 | va_list args;
|
---|
1410 |
|
---|
1411 | // Display the name of the function that called error
|
---|
1412 | printf("[%s]: ", funct);
|
---|
1413 |
|
---|
1414 | // Display the remainder of the message
|
---|
1415 | va_start(args, fmt);
|
---|
1416 | vprintf(fmt, args);
|
---|
1417 | va_end(args);
|
---|
1418 | }
|
---|
1419 | //!@}
|
---|
1420 |
|
---|
1421 |
|
---|
1422 | //!-----------------------------------------------------------
|
---|
1423 | // @name error
|
---|
1424 | //
|
---|
1425 | // @desc function to send an error message, and abort the program
|
---|
1426 | //
|
---|
1427 | // @var funct Name of the caller function
|
---|
1428 | // @var fmt Format to be used (message)
|
---|
1429 | // @var ... Other information to be shown
|
---|
1430 | //
|
---|
1431 | // @date Sat Jun 27 05:58:56 MET DST 1998
|
---|
1432 | //------------------------------------------------------------
|
---|
1433 | // @function
|
---|
1434 |
|
---|
1435 | //!@{
|
---|
1436 | void
|
---|
1437 | error(const char *funct, char *fmt, ...)
|
---|
1438 | {
|
---|
1439 | va_list args;
|
---|
1440 |
|
---|
1441 | // Display the name of the function that called error
|
---|
1442 | fprintf(stdout, "ERROR in %s: ", funct);
|
---|
1443 |
|
---|
1444 | // Display the remainder of the message
|
---|
1445 | va_start(args, fmt);
|
---|
1446 | vfprintf(stdout, fmt, args);
|
---|
1447 | va_end(args);
|
---|
1448 |
|
---|
1449 | perror(funct);
|
---|
1450 |
|
---|
1451 | exit(1);
|
---|
1452 | }
|
---|
1453 | //!@}
|
---|
1454 |
|
---|
1455 |
|
---|
1456 | //!-----------------------------------------------------------
|
---|
1457 | // @name isA
|
---|
1458 | //
|
---|
1459 | // @desc returns TRUE(FALSE), if the flag is(is not) the given
|
---|
1460 | //
|
---|
1461 | // @var s1 String to be searched
|
---|
1462 | // @var flag Flag to compare with string s1
|
---|
1463 | // @return TRUE: both strings match; FALSE: oth.
|
---|
1464 | //
|
---|
1465 | // @date Wed Jul 8 15:25:39 MET DST 1998
|
---|
1466 | //------------------------------------------------------------
|
---|
1467 | // @function
|
---|
1468 |
|
---|
1469 | //!@{
|
---|
1470 | int
|
---|
1471 | isA( char * s1, const char * flag ) {
|
---|
1472 | return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
|
---|
1473 | }
|
---|
1474 | //!@}
|
---|
1475 |
|
---|
1476 |
|
---|
1477 | //!-----------------------------------------------------------
|
---|
1478 | // @name read_ct_file
|
---|
1479 | //
|
---|
1480 | // @desc read CT definition file
|
---|
1481 | //
|
---|
1482 | // @date Sat Jun 27 05:58:56 MET DST 1998
|
---|
1483 | //------------------------------------------------------------
|
---|
1484 | // @function
|
---|
1485 |
|
---|
1486 | //!@{
|
---|
1487 | void
|
---|
1488 | read_ct_file(void)
|
---|
1489 | {
|
---|
1490 | char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
|
---|
1491 | char token[ITEM_MAX_LENGTH]; //@< a single token
|
---|
1492 | int i, j; //@< dummy counters
|
---|
1493 |
|
---|
1494 | log( "read_ct_file", "start.\n" );
|
---|
1495 |
|
---|
1496 | ifstream ctin ( ct_filename );
|
---|
1497 |
|
---|
1498 | if ( ctin.bad() )
|
---|
1499 | error( "read_ct_file",
|
---|
1500 | "Cannot open CT def. file: %s\n", ct_filename );
|
---|
1501 |
|
---|
1502 | // loop till the "end" directive is reached
|
---|
1503 |
|
---|
1504 | while (!ctin.eof()) {
|
---|
1505 |
|
---|
1506 | // get line from stdin
|
---|
1507 |
|
---|
1508 | ctin.getline(line, LINE_MAX_LENGTH);
|
---|
1509 |
|
---|
1510 | // look for each item at the beginning of the line
|
---|
1511 |
|
---|
1512 | for (i=0; i<=define_mirrors; i++)
|
---|
1513 | if (strstr(line, CT_ITEM_NAMES[i]) == line)
|
---|
1514 | break;
|
---|
1515 |
|
---|
1516 | // if it is not a valid line, just ignore it
|
---|
1517 |
|
---|
1518 | if (i == define_mirrors+1)
|
---|
1519 | continue;
|
---|
1520 |
|
---|
1521 | // case block for each directive
|
---|
1522 |
|
---|
1523 | switch ( i ) {
|
---|
1524 |
|
---|
1525 | case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
|
---|
1526 |
|
---|
1527 | // get focal distance
|
---|
1528 |
|
---|
1529 | sscanf(line, "%s %d", token, &ct_Type);
|
---|
1530 |
|
---|
1531 | log( "read_ct_file", "<Type of Telescope>: %s\n",
|
---|
1532 | ((ct_Type==0) ? "CT1" : "MAGIC") );
|
---|
1533 |
|
---|
1534 | break;
|
---|
1535 |
|
---|
1536 | case focal_distance: // <focal distance> [cm]
|
---|
1537 |
|
---|
1538 | // get focal distance
|
---|
1539 |
|
---|
1540 | sscanf(line, "%s %f", token, &ct_Focal_mean);
|
---|
1541 |
|
---|
1542 | log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
|
---|
1543 |
|
---|
1544 | break;
|
---|
1545 |
|
---|
1546 | case focal_std: // s(focal distance) [cm]
|
---|
1547 |
|
---|
1548 | // get focal distance
|
---|
1549 |
|
---|
1550 | sscanf(line, "%s %f", token, &ct_Focal_std);
|
---|
1551 |
|
---|
1552 | log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
|
---|
1553 |
|
---|
1554 | break;
|
---|
1555 |
|
---|
1556 | case point_spread: // <point spread> [cm]
|
---|
1557 |
|
---|
1558 | // get point spread
|
---|
1559 |
|
---|
1560 | sscanf(line, "%s %f", token, &ct_PSpread_mean);
|
---|
1561 |
|
---|
1562 | log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
|
---|
1563 |
|
---|
1564 | break;
|
---|
1565 |
|
---|
1566 | case point_std: // s(point spread) [cm]
|
---|
1567 |
|
---|
1568 | // get point spread
|
---|
1569 |
|
---|
1570 | sscanf(line, "%s %f", token, &ct_PSpread_std);
|
---|
1571 |
|
---|
1572 | log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
|
---|
1573 |
|
---|
1574 | break;
|
---|
1575 |
|
---|
1576 | case adjustment_dev: // s(adjustment_dev) [cm]
|
---|
1577 |
|
---|
1578 | // get point spread
|
---|
1579 |
|
---|
1580 | sscanf(line, "%s %f", token, &ct_Adjustment_std);
|
---|
1581 |
|
---|
1582 | log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
|
---|
1583 |
|
---|
1584 | break;
|
---|
1585 |
|
---|
1586 | case black_spot: // radius of the black spot in the center [cm]
|
---|
1587 |
|
---|
1588 | // get black spot radius
|
---|
1589 |
|
---|
1590 | sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
|
---|
1591 |
|
---|
1592 | log( "read_ct_file", "Radius of the black spots: %f cm\n",
|
---|
1593 | ct_BlackSpot_rad);
|
---|
1594 |
|
---|
1595 | break;
|
---|
1596 |
|
---|
1597 | case r_mirror: // radius of the mirrors [cm]
|
---|
1598 |
|
---|
1599 | // get radius of mirror
|
---|
1600 |
|
---|
1601 | sscanf(line, "%s %f", token, &ct_RMirror);
|
---|
1602 |
|
---|
1603 | log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
|
---|
1604 |
|
---|
1605 | break;
|
---|
1606 |
|
---|
1607 | case n_mirrors: // number of mirrors
|
---|
1608 |
|
---|
1609 | // get the name of the output_file from the line
|
---|
1610 |
|
---|
1611 | sscanf(line, "%s %d", token, &ct_NMirrors);
|
---|
1612 |
|
---|
1613 | log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
|
---|
1614 |
|
---|
1615 | break;
|
---|
1616 |
|
---|
1617 | case camera_width: // camera width [cm]
|
---|
1618 |
|
---|
1619 | // get the name of the ct_file from the line
|
---|
1620 |
|
---|
1621 | sscanf(line, "%s %f", token, &ct_CameraWidth);
|
---|
1622 |
|
---|
1623 | log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
|
---|
1624 |
|
---|
1625 | break;
|
---|
1626 |
|
---|
1627 | case n_pixels: // number of pixels
|
---|
1628 |
|
---|
1629 | // get the name of the output_file from the line
|
---|
1630 |
|
---|
1631 | sscanf(line, "%s %d", token, &ct_NPixels);
|
---|
1632 |
|
---|
1633 | log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
|
---|
1634 |
|
---|
1635 | break;
|
---|
1636 |
|
---|
1637 | case n_centralpixels: // number of central pixels
|
---|
1638 |
|
---|
1639 | // get the name of the output_file from the line
|
---|
1640 |
|
---|
1641 | sscanf(line, "%s %d", token, &ct_NCentralPixels);
|
---|
1642 |
|
---|
1643 | log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
|
---|
1644 |
|
---|
1645 | break;
|
---|
1646 |
|
---|
1647 | case n_gappixels: // number of gap pixels
|
---|
1648 |
|
---|
1649 | // get the name of the output_file from the line
|
---|
1650 |
|
---|
1651 | sscanf(line, "%s %d", token, &ct_NGapPixels);
|
---|
1652 |
|
---|
1653 | log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
|
---|
1654 |
|
---|
1655 | break;
|
---|
1656 |
|
---|
1657 | case pixel_width: // pixel width [cm]
|
---|
1658 |
|
---|
1659 | // get the name of the ct_file from the line
|
---|
1660 |
|
---|
1661 | sscanf(line, "%s %f", token, &ct_PixelWidth);
|
---|
1662 |
|
---|
1663 | ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
|
---|
1664 | ct_PixelWidth_corner_2_corner_half =
|
---|
1665 | ct_PixelWidth_corner_2_corner * 0.50;
|
---|
1666 | ct_Apot = ct_PixelWidth / 2;
|
---|
1667 | ct_2Apot = ct_Apot * 2.0;
|
---|
1668 |
|
---|
1669 | log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
|
---|
1670 |
|
---|
1671 | break;
|
---|
1672 |
|
---|
1673 | case define_mirrors: // read table with the parameters of the mirrors
|
---|
1674 |
|
---|
1675 | log( "read_ct_file", "Table of mirrors data:\n" );
|
---|
1676 |
|
---|
1677 | // check whether the number of mirrors was already set
|
---|
1678 |
|
---|
1679 | if ( ct_NMirrors == 0 )
|
---|
1680 | error( "read_ct_file", "NMirrors was not set.\n" );
|
---|
1681 |
|
---|
1682 | // allocate memory for paths list
|
---|
1683 |
|
---|
1684 | log( "read_ct_file", "Allocating memory for ct_data\n" );
|
---|
1685 |
|
---|
1686 | ct_data = new float*[ct_NMirrors];
|
---|
1687 |
|
---|
1688 | for (i=0; i<ct_NMirrors; i++)
|
---|
1689 | ct_data[i] = new float[CT_NDATA];
|
---|
1690 |
|
---|
1691 | // read data
|
---|
1692 |
|
---|
1693 | log( "read_ct_file", "Reading mirrors data...\n" );
|
---|
1694 |
|
---|
1695 | for (i=0; i<ct_NMirrors; i++)
|
---|
1696 | for (j=0; j<CT_NDATA; j++)
|
---|
1697 | ctin >> ct_data[i][j];
|
---|
1698 |
|
---|
1699 | break;
|
---|
1700 |
|
---|
1701 | } // switch ( i )
|
---|
1702 |
|
---|
1703 | } // end while
|
---|
1704 |
|
---|
1705 | // end
|
---|
1706 |
|
---|
1707 | log( "read_ct_file", "done.\n" );
|
---|
1708 |
|
---|
1709 | return;
|
---|
1710 | }
|
---|
1711 | //!@}
|
---|
1712 |
|
---|
1713 |
|
---|
1714 | //!-----------------------------------------------------------
|
---|
1715 | // @name read_pixels
|
---|
1716 | //
|
---|
1717 | // @desc read pixels data
|
---|
1718 | //
|
---|
1719 | // @date Fri Mar 12 16:33:34 MET 1999
|
---|
1720 | //------------------------------------------------------------
|
---|
1721 | // @function
|
---|
1722 |
|
---|
1723 | //!@{
|
---|
1724 | void
|
---|
1725 | read_pixels(struct camera *pcam)
|
---|
1726 | {
|
---|
1727 | ifstream qefile;
|
---|
1728 | char line[LINE_MAX_LENGTH];
|
---|
1729 | int n, i, j, icount;
|
---|
1730 | float qe;
|
---|
1731 |
|
---|
1732 | //------------------------------------------------------------
|
---|
1733 | // first, pixels' coordinates
|
---|
1734 |
|
---|
1735 | pcam->inumpixels = ct_NPixels;
|
---|
1736 | pcam->inumcentralpixels = ct_NCentralPixels;
|
---|
1737 | pcam->inumgappixels = ct_NGapPixels;
|
---|
1738 | pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
|
---|
1739 | pcam->dpixdiameter_cm = ct_PixelWidth;
|
---|
1740 |
|
---|
1741 | // initialize pixel numbers
|
---|
1742 |
|
---|
1743 | pixary = new float* [2*ct_NCentralPixels];
|
---|
1744 | for ( i=0; i<2*ct_NCentralPixels; ++i )
|
---|
1745 | pixary[i] = new float[2];
|
---|
1746 |
|
---|
1747 | pixneig = new int* [ct_NCentralPixels];
|
---|
1748 | for ( i=0; i<ct_NCentralPixels; ++i ) {
|
---|
1749 | pixneig[i] = new int[6];
|
---|
1750 | for ( j=0; j<6; ++j )
|
---|
1751 | pixneig[i][j] = -1;
|
---|
1752 | }
|
---|
1753 |
|
---|
1754 | npixneig = new int[ct_NCentralPixels];
|
---|
1755 | for ( i=0; i<ct_NCentralPixels; ++i )
|
---|
1756 | npixneig[i] = 0;
|
---|
1757 |
|
---|
1758 | // generate all coordinates
|
---|
1759 |
|
---|
1760 | igen_pixel_coordinates(pcam);
|
---|
1761 |
|
---|
1762 |
|
---|
1763 | // calculate tables of neighbours
|
---|
1764 |
|
---|
1765 | #ifdef __DEBUG__
|
---|
1766 | for ( n=0 ; n<ct_NPixels ; ++n ) {
|
---|
1767 | cout << "Para el pixel " << n << ": ";
|
---|
1768 | for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
|
---|
1769 | if ( pixels_are_neig(n,i) == TRUE ) {
|
---|
1770 | pixneig[n][npixneig[n]] = i;
|
---|
1771 | pixneig[i][npixneig[i]] = n;
|
---|
1772 | cout << i << ' ';
|
---|
1773 | ++npixneig[n];
|
---|
1774 | ++npixneig[i];
|
---|
1775 | }
|
---|
1776 | }
|
---|
1777 | cout << endl << flush;
|
---|
1778 | }
|
---|
1779 | #else // ! __DEBUG__
|
---|
1780 | for ( n=0 ; n<ct_NCentralPixels ; ++n )
|
---|
1781 | for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
|
---|
1782 | if ( pixels_are_neig(n,i) == TRUE ) {
|
---|
1783 | pixneig[n][npixneig[n]] = i;
|
---|
1784 | pixneig[i][npixneig[i]] = n;
|
---|
1785 | ++npixneig[n];
|
---|
1786 | ++npixneig[i];
|
---|
1787 | }
|
---|
1788 | #endif // ! __DEBUG__
|
---|
1789 |
|
---|
1790 | #ifdef __DEBUG__
|
---|
1791 | for ( n=0 ; n<ct_NPixels ; ++n ) {
|
---|
1792 | cout << n << ':';
|
---|
1793 | for ( j=0; j<npixneig[n]; ++j)
|
---|
1794 | cout << ' ' << pixneig[n][j];
|
---|
1795 | cout << endl << flush;
|
---|
1796 | }
|
---|
1797 | #endif // __DEBUG__
|
---|
1798 |
|
---|
1799 | //------------------------------------------------------------
|
---|
1800 | // second, pixels' QE
|
---|
1801 |
|
---|
1802 | // try to open the file
|
---|
1803 |
|
---|
1804 | log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
|
---|
1805 |
|
---|
1806 | qefile.open( QE_FILE );
|
---|
1807 |
|
---|
1808 | // if it is wrong or does not exist, exit
|
---|
1809 |
|
---|
1810 | if ( qefile.bad() )
|
---|
1811 | error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
|
---|
1812 |
|
---|
1813 | // read file
|
---|
1814 |
|
---|
1815 | log("read_pixels", "Reading data . . .\n");
|
---|
1816 |
|
---|
1817 | i=-1;
|
---|
1818 | icount = 0;
|
---|
1819 |
|
---|
1820 | while ( ! qefile.eof() ) {
|
---|
1821 |
|
---|
1822 | // get line from the file
|
---|
1823 |
|
---|
1824 | qefile.getline(line, LINE_MAX_LENGTH);
|
---|
1825 |
|
---|
1826 | // skip if comment
|
---|
1827 |
|
---|
1828 | if ( *line == '#' )
|
---|
1829 | continue;
|
---|
1830 |
|
---|
1831 | // if it is the first valid value, it is the number of QE data points
|
---|
1832 |
|
---|
1833 | if ( i < 0 ) {
|
---|
1834 |
|
---|
1835 | // get the number of datapoints
|
---|
1836 |
|
---|
1837 | sscanf(line, "%d", &pointsQE);
|
---|
1838 |
|
---|
1839 | // allocate memory for the table of QEs
|
---|
1840 |
|
---|
1841 | QE = new float ** [ct_NPixels];
|
---|
1842 |
|
---|
1843 | for ( i=0; i<ct_NPixels; ++i ) {
|
---|
1844 | QE[i] = new float * [2];
|
---|
1845 | QE[i][0] = new float[pointsQE];
|
---|
1846 | QE[i][1] = new float[pointsQE];
|
---|
1847 | }
|
---|
1848 |
|
---|
1849 | QElambda = new float [pointsQE];
|
---|
1850 |
|
---|
1851 | for ( i=0; i<pointsQE; ++i ) {
|
---|
1852 | qefile.getline(line, LINE_MAX_LENGTH);
|
---|
1853 | sscanf(line, "%f", &QElambda[i]);
|
---|
1854 | }
|
---|
1855 |
|
---|
1856 | i=0;
|
---|
1857 |
|
---|
1858 | continue;
|
---|
1859 | }
|
---|
1860 |
|
---|
1861 | // get the values (num-pixel, num-datapoint, QE-value)
|
---|
1862 |
|
---|
1863 | if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
|
---|
1864 | break;
|
---|
1865 |
|
---|
1866 | if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
|
---|
1867 | ((j-1) < pointsQE) && ((j-1) > -1) ) {
|
---|
1868 | QE[i-1][0][j-1] = QElambda[j-1];
|
---|
1869 | QE[i-1][1][j-1] = qe;
|
---|
1870 | }
|
---|
1871 |
|
---|
1872 | if ( i > ct_NPixels) break;
|
---|
1873 |
|
---|
1874 | icount++;
|
---|
1875 |
|
---|
1876 | }
|
---|
1877 |
|
---|
1878 | if(icount/pointsQE < ct_NPixels){
|
---|
1879 | error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",
|
---|
1880 | icount/pointsQE, ct_NPixels );
|
---|
1881 | }
|
---|
1882 |
|
---|
1883 | // close file
|
---|
1884 |
|
---|
1885 | qefile.close();
|
---|
1886 |
|
---|
1887 | // test QE
|
---|
1888 |
|
---|
1889 | for(icount=0; icount< ct_NPixels; icount++){
|
---|
1890 | for(i=0; i<pointsQE; i++){
|
---|
1891 | if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
|
---|
1892 | QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
|
---|
1893 | error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",
|
---|
1894 | icount, i, QE[icount][0][i], QE[icount][1][i] );
|
---|
1895 | }
|
---|
1896 | }
|
---|
1897 | }
|
---|
1898 |
|
---|
1899 | // end
|
---|
1900 |
|
---|
1901 | log("read_pixels", "Done.\n");
|
---|
1902 |
|
---|
1903 | }
|
---|
1904 | //!@}
|
---|
1905 |
|
---|
1906 |
|
---|
1907 | //!-----------------------------------------------------------
|
---|
1908 | // @name pixels_are_neig
|
---|
1909 | //
|
---|
1910 | // @desc check whether two pixels are neighbours
|
---|
1911 | //
|
---|
1912 | // @var pix1 Number of the first pixel
|
---|
1913 | // @var pix2 Number of the second pixel
|
---|
1914 | // @return TRUE: both pixels are neighbours; FALSE: oth.
|
---|
1915 | //
|
---|
1916 | // @date Wed Sep 9 17:58:37 MET DST 1998
|
---|
1917 | //------------------------------------------------------------
|
---|
1918 | // @function
|
---|
1919 |
|
---|
1920 | //!@{
|
---|
1921 | int
|
---|
1922 | pixels_are_neig(int pix1, int pix2)
|
---|
1923 | {
|
---|
1924 | if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
|
---|
1925 | SQR( pixary[pix1][1] - pixary[pix2][1] ) )
|
---|
1926 | > ct_PixelWidth_corner_2_corner )
|
---|
1927 | return ( FALSE );
|
---|
1928 | else
|
---|
1929 | return ( TRUE );
|
---|
1930 | }
|
---|
1931 | //!@}
|
---|
1932 |
|
---|
1933 | //!-----------------------------------------------------------
|
---|
1934 | // @name igen_pixel_coordinates
|
---|
1935 | //
|
---|
1936 | // @desc generate the pixel center coordinates
|
---|
1937 | //
|
---|
1938 | // @var *pcam structure camera containing all the
|
---|
1939 | // camera information
|
---|
1940 | // @return total number of pixels
|
---|
1941 | //
|
---|
1942 | // DP
|
---|
1943 | //
|
---|
1944 | // @date Thu Oct 14 10:41:03 CEST 1999
|
---|
1945 | //------------------------------------------------------------
|
---|
1946 | // @function
|
---|
1947 |
|
---|
1948 | //!@{
|
---|
1949 | /******** igen_pixel_coordinates() *********************************/
|
---|
1950 |
|
---|
1951 | int igen_pixel_coordinates(struct camera *pcam) {
|
---|
1952 | /* generate pixel coordinates, return value is number of pixels */
|
---|
1953 |
|
---|
1954 | int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
|
---|
1955 | float fsegment_fract;
|
---|
1956 | double dtsize;
|
---|
1957 | double dhsize;
|
---|
1958 | double dpsize;
|
---|
1959 | double dxfirst_pix;
|
---|
1960 | double dyfirst_pix;
|
---|
1961 | double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
|
---|
1962 | double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
|
---|
1963 |
|
---|
1964 |
|
---|
1965 | double dstartx, dstarty; /* for the gap pixels and outer pixels */
|
---|
1966 | int j, nrow;
|
---|
1967 |
|
---|
1968 | dpsize = pcam->dpixdiameter_cm;
|
---|
1969 | dtsize = dpsize * sqrt(3.) / 2.;
|
---|
1970 | dhsize = dpsize / 2.;
|
---|
1971 |
|
---|
1972 | /* Loop over central pixels to generate co-ordinates */
|
---|
1973 |
|
---|
1974 | for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
|
---|
1975 |
|
---|
1976 | /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
|
---|
1977 |
|
---|
1978 | pcam->dpixsizefactor[ipixno-1] = 1.;
|
---|
1979 |
|
---|
1980 | in = 0;
|
---|
1981 |
|
---|
1982 | i = 0;
|
---|
1983 | itot_inside_ring = 0;
|
---|
1984 | iring_no = 0;
|
---|
1985 |
|
---|
1986 | /* Calculate the number of pixels out to and including the ring containing pixel number */
|
---|
1987 | /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
|
---|
1988 |
|
---|
1989 | while (itot_inside_ring == 0){
|
---|
1990 |
|
---|
1991 | iN = 3*(i*(i+1)) + 1;
|
---|
1992 |
|
---|
1993 | if (ipixno <= iN){
|
---|
1994 | iring_no = i;
|
---|
1995 | itot_inside_ring = iN;
|
---|
1996 | }
|
---|
1997 |
|
---|
1998 | i++;
|
---|
1999 | }
|
---|
2000 |
|
---|
2001 |
|
---|
2002 | /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
|
---|
2003 |
|
---|
2004 | ipix_in_ring = 0;
|
---|
2005 | for (i = 0; i < iring_no; ++i){
|
---|
2006 |
|
---|
2007 | ipix_in_ring = ipix_in_ring + 6;
|
---|
2008 | }
|
---|
2009 |
|
---|
2010 | /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
|
---|
2011 | /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
|
---|
2012 | /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
|
---|
2013 |
|
---|
2014 | isegment = 0;
|
---|
2015 | fsegment_fract = 0.;
|
---|
2016 | if (iring_no > 0) {
|
---|
2017 |
|
---|
2018 | isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
|
---|
2019 |
|
---|
2020 | fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
|
---|
2021 |
|
---|
2022 | }
|
---|
2023 |
|
---|
2024 | /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
|
---|
2025 | /* pixel width (flat to flat) dpsize. */
|
---|
2026 |
|
---|
2027 | dxfirst_pix = dpsize*iring_no;
|
---|
2028 | dyfirst_pix = 0.;
|
---|
2029 |
|
---|
2030 | /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
|
---|
2031 |
|
---|
2032 | ddxseg1 = - dhsize*iring_no;
|
---|
2033 | ddyseg1 = dtsize*iring_no;
|
---|
2034 | ddxseg2 = -dpsize*iring_no;
|
---|
2035 | ddyseg2 = 0.;
|
---|
2036 | ddxseg3 = ddxseg1;
|
---|
2037 | ddyseg3 = -ddyseg1;
|
---|
2038 | ddxseg4 = -ddxseg1;
|
---|
2039 | ddyseg4 = -ddyseg1;
|
---|
2040 | ddxseg5 = -ddxseg2;
|
---|
2041 | ddyseg5 = 0.;
|
---|
2042 | ddxseg6 = -ddxseg1;
|
---|
2043 | ddyseg6 = ddyseg1;
|
---|
2044 |
|
---|
2045 | /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
|
---|
2046 | /* anti-clockwise around the ring by adding the segment to segment vectors. */
|
---|
2047 |
|
---|
2048 | switch (isegment) {
|
---|
2049 |
|
---|
2050 | case 0:
|
---|
2051 |
|
---|
2052 | pcam->dxc[ipixno-1] = 0.;
|
---|
2053 | pcam->dyc[ipixno-1] = 0.;
|
---|
2054 |
|
---|
2055 | case 1:
|
---|
2056 | pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
|
---|
2057 | pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
|
---|
2058 |
|
---|
2059 | break;
|
---|
2060 |
|
---|
2061 | case 2:
|
---|
2062 |
|
---|
2063 | pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
|
---|
2064 | pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
|
---|
2065 |
|
---|
2066 | break;
|
---|
2067 |
|
---|
2068 | case 3:
|
---|
2069 |
|
---|
2070 | pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
|
---|
2071 | pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
|
---|
2072 |
|
---|
2073 | break;
|
---|
2074 |
|
---|
2075 | case 4:
|
---|
2076 |
|
---|
2077 | pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
|
---|
2078 | pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
|
---|
2079 |
|
---|
2080 | break;
|
---|
2081 |
|
---|
2082 | case 5:
|
---|
2083 |
|
---|
2084 | pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
|
---|
2085 | pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
|
---|
2086 |
|
---|
2087 | break;
|
---|
2088 |
|
---|
2089 | case 6:
|
---|
2090 |
|
---|
2091 | pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
|
---|
2092 | pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
|
---|
2093 |
|
---|
2094 | break;
|
---|
2095 |
|
---|
2096 | default:
|
---|
2097 |
|
---|
2098 | fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
|
---|
2099 | return(0);
|
---|
2100 |
|
---|
2101 | } /* end switch */
|
---|
2102 |
|
---|
2103 | } /* end for */
|
---|
2104 |
|
---|
2105 | dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
|
---|
2106 | dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
|
---|
2107 |
|
---|
2108 | if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
|
---|
2109 |
|
---|
2110 | j = pcam->inumcentralpixels;
|
---|
2111 |
|
---|
2112 | for(i=0; i<pcam->inumgappixels; i=i+6){
|
---|
2113 | pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
|
---|
2114 | pcam->dyc[j + i ] = dstarty;
|
---|
2115 | pcam->dpixsizefactor[j + i] = 1.;
|
---|
2116 | pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
|
---|
2117 | pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
|
---|
2118 | pcam->dpixsizefactor[j + i + 1] = 1.;
|
---|
2119 | pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
|
---|
2120 | pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
|
---|
2121 | pcam->dpixsizefactor[j + i+ 2] = 1.;
|
---|
2122 | pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
|
---|
2123 | pcam->dyc[j + i + 3] = dstarty;
|
---|
2124 | pcam->dpixsizefactor[j + i+ 3] = 1.;
|
---|
2125 | pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
|
---|
2126 | pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
|
---|
2127 | pcam->dpixsizefactor[j + i+ 4] = 1.;
|
---|
2128 | pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
|
---|
2129 | pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
|
---|
2130 | pcam->dpixsizefactor[j + i + 5] = 1.;
|
---|
2131 | } /* end for */
|
---|
2132 | } /* end if */
|
---|
2133 |
|
---|
2134 | /* generate positions of the outer pixels */
|
---|
2135 |
|
---|
2136 | if( pcam->inumbigpixels > 0 ){
|
---|
2137 |
|
---|
2138 | j = pcam->inumcentralpixels + pcam->inumgappixels;
|
---|
2139 |
|
---|
2140 | for(i=0; i<pcam->inumbigpixels; i++){
|
---|
2141 | pcam->dpixsizefactor[j + i] = 2.;
|
---|
2142 | }
|
---|
2143 |
|
---|
2144 | in = 0;
|
---|
2145 |
|
---|
2146 | nrow = (int) ceil(dstartx / 2. / dpsize);
|
---|
2147 |
|
---|
2148 | while(in < pcam->inumbigpixels){
|
---|
2149 |
|
---|
2150 | pcam->dxc[j + in] = dstartx + dpsize;
|
---|
2151 | pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
|
---|
2152 | pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
|
---|
2153 | pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
|
---|
2154 | pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
|
---|
2155 | pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
|
---|
2156 | pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
|
---|
2157 | pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
|
---|
2158 | pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
|
---|
2159 | pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
|
---|
2160 | pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
|
---|
2161 | pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
|
---|
2162 | for(i=1; i<nrow; i++){
|
---|
2163 | pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
|
---|
2164 | pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
|
---|
2165 | pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
|
---|
2166 | pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
|
---|
2167 | pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
|
---|
2168 | pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
|
---|
2169 | pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
|
---|
2170 | pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
|
---|
2171 | pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
|
---|
2172 | pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
|
---|
2173 | pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
|
---|
2174 | pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
|
---|
2175 | }
|
---|
2176 | in = in + 6 * nrow;
|
---|
2177 | dstartx = dstartx + 2. * dpsize;
|
---|
2178 | nrow = nrow + 1;
|
---|
2179 | } /* end while */
|
---|
2180 |
|
---|
2181 | } /* end if */
|
---|
2182 |
|
---|
2183 | return(pcam->inumpixels);
|
---|
2184 |
|
---|
2185 | }
|
---|
2186 | //!@}
|
---|
2187 |
|
---|
2188 | //!-----------------------------------------------------------
|
---|
2189 | // @name bpoint_is_in_pix
|
---|
2190 | //
|
---|
2191 | // @desc check if a point (x,y) in camera coordinates is inside a given pixel
|
---|
2192 | //
|
---|
2193 | // @var *pcam structure camera containing all the
|
---|
2194 | // camera information
|
---|
2195 | // @var dx, dy point coordinates in centimeters
|
---|
2196 | // @var ipixnum pixel number (starting at 0)
|
---|
2197 | // @return TRUE if the point is inside the pixel, FALSE otherwise
|
---|
2198 | //
|
---|
2199 | // DP
|
---|
2200 | //
|
---|
2201 | // @date Thu Oct 14 16:59:04 CEST 1999
|
---|
2202 | //------------------------------------------------------------
|
---|
2203 | // @function
|
---|
2204 |
|
---|
2205 | //!@{
|
---|
2206 |
|
---|
2207 | /******** bpoint_is_in_pix() ***************************************/
|
---|
2208 |
|
---|
2209 | int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
|
---|
2210 | /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
|
---|
2211 | /* the pixel is assumed to be a "closed set" */
|
---|
2212 |
|
---|
2213 | double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
|
---|
2214 | double c, xx, yy; /* auxiliary variable */
|
---|
2215 |
|
---|
2216 | b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
|
---|
2217 | a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
|
---|
2218 | c = 1. - 1./sqrt(3.);
|
---|
2219 | if((ipixnum < 0)||(ipixnum >= pcam->inumpixels)){
|
---|
2220 | fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
|
---|
2221 | fprintf(stderr, "Exiting.\n");
|
---|
2222 | exit(203);
|
---|
2223 | }
|
---|
2224 | xx = dx - pcam->dxc[ipixnum];
|
---|
2225 | yy = dy - pcam->dyc[ipixnum];
|
---|
2226 |
|
---|
2227 | if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) ||
|
---|
2228 | ((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
|
---|
2229 | return(TRUE); /* point is inside */
|
---|
2230 | }
|
---|
2231 | else{
|
---|
2232 | return(FALSE); /* point is outside */
|
---|
2233 | }
|
---|
2234 | }
|
---|
2235 |
|
---|
2236 | //!@}
|
---|
2237 |
|
---|
2238 | //------------------------------------------------------------
|
---|
2239 | // @name dist_r_P
|
---|
2240 | //
|
---|
2241 | // @desc distance straight line r - point P
|
---|
2242 | //
|
---|
2243 | // @date Sat Jun 27 05:58:56 MET DST 1998
|
---|
2244 | // @function @code
|
---|
2245 | //------------------------------------------------------------
|
---|
2246 | // dist_r_P
|
---|
2247 | //
|
---|
2248 | // distance straight line r - point P
|
---|
2249 | //------------------------------------------------------------
|
---|
2250 |
|
---|
2251 | float
|
---|
2252 | dist_r_P(float a, float b, float c,
|
---|
2253 | float l, float m, float n,
|
---|
2254 | float x, float y, float z)
|
---|
2255 | {
|
---|
2256 | return (
|
---|
2257 | sqrt((SQR((a-x)*m-(b-y)*l) +
|
---|
2258 | SQR((b-y)*n-(c-z)*m) +
|
---|
2259 | SQR((c-z)*l-(a-x)*n))/
|
---|
2260 | (SQR(l)+SQR(m)+SQR(n))
|
---|
2261 | )
|
---|
2262 | );
|
---|
2263 | }
|
---|
2264 |
|
---|
2265 | //------------------------------------------------------------
|
---|
2266 | // @name check_reflector_file
|
---|
2267 | //
|
---|
2268 | // @desc check if a given reflector file has the right signature
|
---|
2269 | // @desc return TRUE or FALSE
|
---|
2270 | //
|
---|
2271 | // @date Mon Feb 14 16:44:21 CET 2000
|
---|
2272 | // @function @code
|
---|
2273 | //------------------------------------------------------------
|
---|
2274 |
|
---|
2275 | int check_reflector_file(FILE *infile){
|
---|
2276 |
|
---|
2277 | char Signature[20]; // auxiliary variable
|
---|
2278 | char sign[20]; // auxiliary variable
|
---|
2279 |
|
---|
2280 | strcpy(Signature, REFL_SIGNATURE);
|
---|
2281 |
|
---|
2282 | strcpy(sign, Signature);
|
---|
2283 |
|
---|
2284 | fread( (char *)sign, strlen(Signature), 1, infile);
|
---|
2285 |
|
---|
2286 | if (strcmp(sign, Signature) != 0) {
|
---|
2287 | cout << "ERROR: Signature of .rfl file is not correct\n";
|
---|
2288 | cout << '"' << sign << '"' << '\n';
|
---|
2289 | cout << "should be: " << Signature << '\n';
|
---|
2290 | return(FALSE);
|
---|
2291 | }
|
---|
2292 |
|
---|
2293 | fread( (char *)sign, 1, 1, infile);
|
---|
2294 |
|
---|
2295 | return(TRUE);
|
---|
2296 |
|
---|
2297 | }
|
---|
2298 |
|
---|
2299 | //------------------------------------------------------------
|
---|
2300 | // @name lin_interpol
|
---|
2301 | //
|
---|
2302 | // @desc interpolate linearly between two points returning the
|
---|
2303 | // @desc y-value of the result
|
---|
2304 | //
|
---|
2305 | // @date Thu Feb 17 11:31:32 CET 2000
|
---|
2306 | // @function @code
|
---|
2307 | //------------------------------------------------------------
|
---|
2308 |
|
---|
2309 | float lin_interpol( float x1, float y1, float x2, float y2, float x){
|
---|
2310 |
|
---|
2311 | if( (x2 - x1)<=0. ){ // avoid division by zero, return average
|
---|
2312 | cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
|
---|
2313 | return((y1+y2)/2.);
|
---|
2314 | }
|
---|
2315 | else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
|
---|
2316 | return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
|
---|
2317 | }
|
---|
2318 | }
|
---|
2319 |
|
---|
2320 |
|
---|
2321 | //------------------------------------------------------------
|
---|
2322 | // @name Photoelectron
|
---|
2323 | //
|
---|
2324 | // @desc constructor for class Photoelectron
|
---|
2325 | //
|
---|
2326 | // @date Mon Feb 15 16:44:21 CET 2000
|
---|
2327 | // @function @code
|
---|
2328 | //------------------------------------------------------------
|
---|
2329 |
|
---|
2330 | Photoelectron::Photoelectron(void){
|
---|
2331 | iarrtime_ns = NOTIME;
|
---|
2332 | ipixnum = -1;
|
---|
2333 | }
|
---|
2334 |
|
---|
2335 | //------------------------------------------------------------
|
---|
2336 | // @name produce_phes
|
---|
2337 | //
|
---|
2338 | // @desc read the photons of an event, pixelize them and simulate QE
|
---|
2339 | // @desc return various statistics and the array of Photoelectrons
|
---|
2340 | //
|
---|
2341 | // @date Mon Feb 14 16:44:21 CET 2000
|
---|
2342 | // @function @code
|
---|
2343 | //------------------------------------------------------------
|
---|
2344 |
|
---|
2345 | int produce_phes( FILE *sp, // the input file
|
---|
2346 | struct camera *cam, // the camera layout
|
---|
2347 | float minwl_nm, // the minimum accepted wavelength
|
---|
2348 | float maxwl_nm, // the maximum accepted wavelength
|
---|
2349 | class Photoelectron phe[iMAXNUMPHE], // the generated phes
|
---|
2350 | int *itotnphe, // total number of produced photoelectrons
|
---|
2351 | float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
|
---|
2352 | int *incph, // total number of cph read
|
---|
2353 | float *tmin_ns, // minimum arrival time of all phes
|
---|
2354 | float *tmax_ns // maximum arrival time of all phes
|
---|
2355 | ){
|
---|
2356 |
|
---|
2357 | static int i, k, ipixnum;
|
---|
2358 | static float cx, cy, wl, qe, t;
|
---|
2359 | static MCCphoton photon;
|
---|
2360 | static float **qept;
|
---|
2361 | static char flag[SIZE_OF_FLAGS + 1];
|
---|
2362 | static float radius;
|
---|
2363 |
|
---|
2364 | // reset variables
|
---|
2365 |
|
---|
2366 | for ( i=0; i<cam->inumpixels; ++i ){
|
---|
2367 |
|
---|
2368 | nphe[i] = 0.0;
|
---|
2369 |
|
---|
2370 | }
|
---|
2371 |
|
---|
2372 | *itotnphe = 0;
|
---|
2373 | *incph = 0;
|
---|
2374 | *tmin_ns = NOTIME; // very big
|
---|
2375 | *tmax_ns = -NOTIME; // very small
|
---|
2376 |
|
---|
2377 | radius = cam->dxc[cam->inumpixels-1]
|
---|
2378 | + 1.5*cam->dpixdiameter_cm*cam->dpixsizefactor[cam->inumpixels-1];
|
---|
2379 |
|
---|
2380 | //- - - - - - - - - - - - - - - - - - - - - - - - -
|
---|
2381 | // read photons and "map" them into the pixels
|
---|
2382 | //--------------------------------------------------
|
---|
2383 |
|
---|
2384 | // initialize CPhoton
|
---|
2385 |
|
---|
2386 | photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
|
---|
2387 |
|
---|
2388 | // read the photons data
|
---|
2389 |
|
---|
2390 | fread ( flag, SIZE_OF_FLAGS, 1, sp );
|
---|
2391 |
|
---|
2392 | // loop over the photons
|
---|
2393 |
|
---|
2394 | while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
|
---|
2395 |
|
---|
2396 | memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
|
---|
2397 |
|
---|
2398 | fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
|
---|
2399 |
|
---|
2400 | // increase number of photons
|
---|
2401 |
|
---|
2402 | (*incph)++;
|
---|
2403 |
|
---|
2404 | //
|
---|
2405 | // Pixelization
|
---|
2406 | //
|
---|
2407 |
|
---|
2408 | cx = photon.get_x();
|
---|
2409 | cy = photon.get_y();
|
---|
2410 |
|
---|
2411 | // get wavelength
|
---|
2412 |
|
---|
2413 | wl = photon.get_wl();
|
---|
2414 |
|
---|
2415 | // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n";
|
---|
2416 |
|
---|
2417 | // check if photon has valid wavelength and is inside outermost camera radius
|
---|
2418 |
|
---|
2419 | if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){
|
---|
2420 |
|
---|
2421 | // cout << " lost first check\n";
|
---|
2422 |
|
---|
2423 | // read next CPhoton
|
---|
2424 |
|
---|
2425 | fread ( flag, SIZE_OF_FLAGS, 1, sp );
|
---|
2426 |
|
---|
2427 | // go to beginning of loop, the photon is lost
|
---|
2428 | continue;
|
---|
2429 |
|
---|
2430 | }
|
---|
2431 |
|
---|
2432 | ipixnum = -1;
|
---|
2433 |
|
---|
2434 | for(i=0; i<cam->inumpixels; i++){
|
---|
2435 | if( bpoint_is_in_pix( cx, cy, i, cam) ){
|
---|
2436 | ipixnum = i;
|
---|
2437 | break;
|
---|
2438 | }
|
---|
2439 | }
|
---|
2440 |
|
---|
2441 | if(ipixnum==-1){// the photon is in none of the pixels
|
---|
2442 |
|
---|
2443 | // cout << " lost pixlization\n";
|
---|
2444 |
|
---|
2445 | // read next CPhoton
|
---|
2446 |
|
---|
2447 | fread ( flag, SIZE_OF_FLAGS, 1, sp );
|
---|
2448 |
|
---|
2449 | // go to beginning of loop, the photon is lost
|
---|
2450 | continue;
|
---|
2451 | }
|
---|
2452 |
|
---|
2453 | //+++
|
---|
2454 | // QE simulation
|
---|
2455 | //---
|
---|
2456 |
|
---|
2457 | // set pointer to the QE table of the relevant pixel
|
---|
2458 |
|
---|
2459 | qept = (float **)QE[ipixnum];
|
---|
2460 |
|
---|
2461 | // check if wl is inside table; outside the table, QE is assumed to be zero
|
---|
2462 |
|
---|
2463 | if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
|
---|
2464 |
|
---|
2465 | // cout << " lost\n";
|
---|
2466 |
|
---|
2467 | // read next Photon
|
---|
2468 |
|
---|
2469 | fread ( flag, SIZE_OF_FLAGS, 1, sp );
|
---|
2470 |
|
---|
2471 | // go to beginning of loop
|
---|
2472 | continue;
|
---|
2473 |
|
---|
2474 | }
|
---|
2475 |
|
---|
2476 | // find data point in the QE table (-> k)
|
---|
2477 |
|
---|
2478 | k = 1; // start at 1 because the condition was already tested for 0
|
---|
2479 | while (k < pointsQE-1 && qept[0][k] < wl){
|
---|
2480 | k++;
|
---|
2481 | }
|
---|
2482 |
|
---|
2483 | // calculate the qe value between 0. and 1.
|
---|
2484 |
|
---|
2485 | qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
|
---|
2486 |
|
---|
2487 | // if random > quantum efficiency, reject it
|
---|
2488 |
|
---|
2489 | if ( RandomNumber > qe ) {
|
---|
2490 |
|
---|
2491 | // cout << " lost\n";
|
---|
2492 |
|
---|
2493 | // read next Photon
|
---|
2494 |
|
---|
2495 | fread ( flag, SIZE_OF_FLAGS, 1, sp );
|
---|
2496 |
|
---|
2497 | // go to beginning of loop
|
---|
2498 | continue;
|
---|
2499 |
|
---|
2500 | }
|
---|
2501 |
|
---|
2502 | //+++
|
---|
2503 | // The photon has produced a photo electron
|
---|
2504 | //---
|
---|
2505 |
|
---|
2506 | // cout << " accepted\n";
|
---|
2507 |
|
---|
2508 | // increment the number of photoelectrons in the relevant pixel
|
---|
2509 |
|
---|
2510 | nphe[ipixnum] += 1.0;
|
---|
2511 |
|
---|
2512 | t = photon.get_t() ;
|
---|
2513 |
|
---|
2514 | // cout << " t " << t;
|
---|
2515 |
|
---|
2516 | // find minimum and maximum arrival time
|
---|
2517 |
|
---|
2518 | if(t < *tmin_ns){
|
---|
2519 | *tmin_ns = t; // memorize time
|
---|
2520 | }
|
---|
2521 | if(t > *tmax_ns){
|
---|
2522 | *tmax_ns = t; // memorize time
|
---|
2523 | }
|
---|
2524 |
|
---|
2525 | // store the new photoelectron
|
---|
2526 |
|
---|
2527 | if(*itotnphe >= iMAXNUMPHE){
|
---|
2528 | cout << "Error: Memory overflow. Event produces more than maximum\n";
|
---|
2529 | cout << " allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
|
---|
2530 | return(1);
|
---|
2531 | }
|
---|
2532 |
|
---|
2533 | phe[*itotnphe].iarrtime_ns = (int)t;
|
---|
2534 | phe[*itotnphe].ipixnum = ipixnum;
|
---|
2535 |
|
---|
2536 | *itotnphe += 1;
|
---|
2537 |
|
---|
2538 | // read next Photon
|
---|
2539 |
|
---|
2540 | fread( flag, SIZE_OF_FLAGS, 1, sp );
|
---|
2541 |
|
---|
2542 | } // end while, i.e. found end of event
|
---|
2543 |
|
---|
2544 | return(0);
|
---|
2545 |
|
---|
2546 | }
|
---|
2547 |
|
---|
2548 |
|
---|
2549 | //------------------------------------------------------------
|
---|
2550 | // @name produce_nsbrates
|
---|
2551 | //
|
---|
2552 | // @desc read the starfield file, call produce_phes on it in,
|
---|
2553 | // @desc different wavebands, calculate the nsbrates
|
---|
2554 | //
|
---|
2555 | // @date Mon Feb 14 16:44:21 CET 2000
|
---|
2556 | // @function @code
|
---|
2557 | //------------------------------------------------------------
|
---|
2558 |
|
---|
2559 | int produce_nsbrates( char *iname, // the starfield input file name
|
---|
2560 | struct camera *cam, // camera layout
|
---|
2561 | class Photoelectron phe[iMAXNUMPHE], // array for photoelectrons
|
---|
2562 | float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
|
---|
2563 | // the NSB rates in phe/ns for each pixel
|
---|
2564 | ){
|
---|
2565 |
|
---|
2566 | int i, j, k, ii; // counters
|
---|
2567 |
|
---|
2568 | static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
|
---|
2569 | WAVEBANDBOUND2,
|
---|
2570 | WAVEBANDBOUND3,
|
---|
2571 | WAVEBANDBOUND4,
|
---|
2572 | WAVEBANDBOUND5,
|
---|
2573 | WAVEBANDBOUND6 };
|
---|
2574 |
|
---|
2575 | FILE *infile; // the input file
|
---|
2576 | fpos_t fileposition; // marker on the input file
|
---|
2577 | static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
|
---|
2578 | static MCEventHeader evth; // the event header
|
---|
2579 | static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel
|
---|
2580 | int itnphe; // total number of produced photoelectrons
|
---|
2581 | int itotnphe; // total number of produced photoelectrons after averaging
|
---|
2582 | int incph; // total number of cph read
|
---|
2583 | float tmin_ns; // minimum arrival time of all phes
|
---|
2584 | float tmax_ns; // maximum arrival time of all phes
|
---|
2585 | float integtime_ns; // integration time from the starfield generator
|
---|
2586 |
|
---|
2587 | // open input file
|
---|
2588 |
|
---|
2589 | log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
|
---|
2590 |
|
---|
2591 | infile = fopen( iname, "r" );
|
---|
2592 | if ( infile == NULL )
|
---|
2593 | error( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
|
---|
2594 |
|
---|
2595 | // get signature, and check it
|
---|
2596 |
|
---|
2597 | if(check_reflector_file(infile)==FALSE){
|
---|
2598 | exit(1);
|
---|
2599 | }
|
---|
2600 |
|
---|
2601 | // initialize flag
|
---|
2602 |
|
---|
2603 | strcpy( cflag, " \0" );
|
---|
2604 |
|
---|
2605 | // get flag
|
---|
2606 |
|
---|
2607 | fread( cflag, SIZE_OF_FLAGS, 1, infile );
|
---|
2608 |
|
---|
2609 | if ( ! feof(infile)){
|
---|
2610 |
|
---|
2611 | // reading .rfl file
|
---|
2612 |
|
---|
2613 | if(!isA( cflag, FLAG_START_OF_RUN )){
|
---|
2614 | error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
|
---|
2615 | }
|
---|
2616 | else { // found start of run
|
---|
2617 |
|
---|
2618 | fread( cflag, SIZE_OF_FLAGS, 1, infile );
|
---|
2619 |
|
---|
2620 | if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event
|
---|
2621 |
|
---|
2622 | // get MCEventHeader
|
---|
2623 |
|
---|
2624 | fread( (char*)&evth, evth.mysize(), 1, infile );
|
---|
2625 |
|
---|
2626 | integtime_ns = evth.get_energy();
|
---|
2627 |
|
---|
2628 | // memorize where we are in the file
|
---|
2629 |
|
---|
2630 | if (fgetpos( infile, &fileposition ) != 0){
|
---|
2631 | error( SIGNATURE, "Cannot position in file ...\n");
|
---|
2632 | }
|
---|
2633 |
|
---|
2634 | // loop over the wavebands
|
---|
2635 |
|
---|
2636 | for(i=0; i<iNUMWAVEBANDS; i++){
|
---|
2637 |
|
---|
2638 | // initialize the rate array
|
---|
2639 |
|
---|
2640 | for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
|
---|
2641 | rate_phepns[j][i] = 0.;
|
---|
2642 | }
|
---|
2643 |
|
---|
2644 | itotnphe = 0;
|
---|
2645 |
|
---|
2646 | // read the photons and produce the photoelectrons
|
---|
2647 | // - in order to average over the QE simulation, call the
|
---|
2648 | // production function iNUMNSBPRODCALLS times
|
---|
2649 |
|
---|
2650 | for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
|
---|
2651 |
|
---|
2652 | // position the file pointer to the beginning of the photons
|
---|
2653 |
|
---|
2654 | fsetpos( infile, &fileposition);
|
---|
2655 |
|
---|
2656 | // produce photoelectrons
|
---|
2657 |
|
---|
2658 | k = produce_phes( infile,
|
---|
2659 | cam,
|
---|
2660 | wl_nm[i],
|
---|
2661 | wl_nm[i+1],
|
---|
2662 | phe, // this is a dummy here
|
---|
2663 | &itnphe,
|
---|
2664 | nphe, // we want this!
|
---|
2665 | &incph,
|
---|
2666 | &tmin_ns,
|
---|
2667 | &tmax_ns );
|
---|
2668 |
|
---|
2669 | if( k != 0 ){ // non-zero returnvalue means error
|
---|
2670 | cout << "Exiting.\n";
|
---|
2671 | exit(1);
|
---|
2672 | }
|
---|
2673 |
|
---|
2674 | for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
|
---|
2675 | rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
|
---|
2676 | }
|
---|
2677 |
|
---|
2678 | itotnphe += itnphe;
|
---|
2679 |
|
---|
2680 | } // end for(ii=0 ...
|
---|
2681 |
|
---|
2682 | fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
|
---|
2683 | wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
|
---|
2684 |
|
---|
2685 | } // end for(i=0 ...
|
---|
2686 |
|
---|
2687 | }
|
---|
2688 | else{
|
---|
2689 | cout << "Starfield file contains no event.\nExiting.\n";
|
---|
2690 | exit(1);
|
---|
2691 | } // end if( isA ... event
|
---|
2692 | } // end if ( isA ... run
|
---|
2693 | }
|
---|
2694 | else{
|
---|
2695 | cout << "Starfield file contains no run.\nExiting.\n";
|
---|
2696 | exit(1);
|
---|
2697 | }
|
---|
2698 |
|
---|
2699 | fclose( infile );
|
---|
2700 |
|
---|
2701 | return(0);
|
---|
2702 |
|
---|
2703 | }
|
---|
2704 |
|
---|
2705 |
|
---|
2706 | //------------------------------------------------------------
|
---|
2707 | // @name produce_nsb_phes
|
---|
2708 | //
|
---|
2709 | // @desc produce the photoelectrons from the nsbrates
|
---|
2710 | //
|
---|
2711 | // @date Thu Feb 17 17:10:40 CET 2000
|
---|
2712 | // @function @code
|
---|
2713 | //------------------------------------------------------------
|
---|
2714 |
|
---|
2715 | int produce_nsb_phes( float *atmin_ns,
|
---|
2716 | float *atmax_ns,
|
---|
2717 | float theta_rad,
|
---|
2718 | struct camera *cam,
|
---|
2719 | float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
|
---|
2720 | float difnsbr_phepns[iMAXNUMPIX],
|
---|
2721 | float extinction[iNUMWAVEBANDS],
|
---|
2722 | float fnpx[iMAXNUMPIX],
|
---|
2723 | Photoelectron photo[iMAXNUMPHE],
|
---|
2724 | int *inphe,
|
---|
2725 | float base_mv[iMAXNUMPIX]){
|
---|
2726 |
|
---|
2727 | float simtime_ns; // NSB simulation time
|
---|
2728 | int i, j, k, ii;
|
---|
2729 | float zenfactor; // correction factor calculated from the extinction
|
---|
2730 | int inumnsbphe; // number of photoelectrons caused by NSB
|
---|
2731 |
|
---|
2732 | ii = *inphe; // avoid dereferencing
|
---|
2733 |
|
---|
2734 | // check if the arrival times are set; if not generate them
|
---|
2735 |
|
---|
2736 | if(*atmin_ns == NOTIME){
|
---|
2737 |
|
---|
2738 | *atmin_ns = 0.;
|
---|
2739 | *atmax_ns = simtime_ns = SLICES*WIDTH_TIMESLICE;
|
---|
2740 |
|
---|
2741 | }
|
---|
2742 | else{ // extend the event time window by the given offsets
|
---|
2743 |
|
---|
2744 | *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
|
---|
2745 | *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
|
---|
2746 |
|
---|
2747 | simtime_ns = *atmax_ns - *atmin_ns;
|
---|
2748 |
|
---|
2749 | // make sure the simulated time is long enough for the FADC simulation
|
---|
2750 |
|
---|
2751 | if(simtime_ns< SLICES*WIDTH_TIMESLICE){
|
---|
2752 | *atmax_ns = *atmin_ns + SLICES*WIDTH_TIMESLICE;
|
---|
2753 | simtime_ns = SLICES*WIDTH_TIMESLICE;
|
---|
2754 | }
|
---|
2755 |
|
---|
2756 | }
|
---|
2757 |
|
---|
2758 | // initialize baselines
|
---|
2759 |
|
---|
2760 | for(i=0; i<cam->inumpixels; i++){
|
---|
2761 | base_mv[i] = 0.;
|
---|
2762 | }
|
---|
2763 |
|
---|
2764 | // calculate baselines and generate phes
|
---|
2765 |
|
---|
2766 | for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
|
---|
2767 |
|
---|
2768 | // calculate the effect of the atmospheric extinction
|
---|
2769 |
|
---|
2770 | zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
|
---|
2771 |
|
---|
2772 | for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
|
---|
2773 |
|
---|
2774 | inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)
|
---|
2775 | * simtime_ns );
|
---|
2776 |
|
---|
2777 | base_mv[j] += inumnsbphe;
|
---|
2778 |
|
---|
2779 | // randomize
|
---|
2780 |
|
---|
2781 | inumnsbphe = ignpoi( inumnsbphe );
|
---|
2782 |
|
---|
2783 | // create the photoelectrons
|
---|
2784 |
|
---|
2785 | for(k=0; k < inumnsbphe; k++){
|
---|
2786 |
|
---|
2787 | if(ii >= iMAXNUMPHE){
|
---|
2788 | cout << "Error: Memory overflow. NSB simulation produces more than maximum\n";
|
---|
2789 | cout << " allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
|
---|
2790 | return(1);
|
---|
2791 | }
|
---|
2792 |
|
---|
2793 | photo[ii].iarrtime_ns = (int)(RandomNumber * simtime_ns + *atmin_ns );
|
---|
2794 | photo[ii].ipixnum = j;
|
---|
2795 |
|
---|
2796 | // cout << "Created phe " << photo[ii].iarrtime_ns << " "
|
---|
2797 | // << photo[ii].ipixnum << "\n";
|
---|
2798 |
|
---|
2799 | ii++; // increment total number of photoelectons
|
---|
2800 |
|
---|
2801 | fnpx[j] += 1.; // increment number of photoelectrons in this pixel
|
---|
2802 |
|
---|
2803 | }
|
---|
2804 |
|
---|
2805 | } // end for(j=0 ...
|
---|
2806 | } // end for(i=0 ...
|
---|
2807 |
|
---|
2808 | // finish baseline calculation
|
---|
2809 |
|
---|
2810 | for(i=0; i<cam->inumpixels; i++){
|
---|
2811 | base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
|
---|
2812 | }
|
---|
2813 |
|
---|
2814 | *inphe = ii; // update the pointer
|
---|
2815 |
|
---|
2816 | return(0);
|
---|
2817 | }
|
---|
2818 |
|
---|
2819 |
|
---|
2820 | // @endcode
|
---|
2821 |
|
---|
2822 |
|
---|
2823 | //=------------------------------------------------------------
|
---|
2824 | //!@subsection Log of this file.
|
---|
2825 |
|
---|
2826 | //!@{
|
---|
2827 | //
|
---|
2828 | // $Log: not supported by cvs2svn $
|
---|
2829 | // Revision 1.5 2000/02/18 17:40:35 petry
|
---|
2830 | // This version includes drastic changes compared to camera.cxx 1.4.
|
---|
2831 | // It is not yet finished and not immediately useful because the
|
---|
2832 | // trigger simulation is not yet re-implemented. I had to take it
|
---|
2833 | // out together with some other stuff in order to tidy the whole
|
---|
2834 | // program up. This is not meant as an insult to anyone. I needed
|
---|
2835 | // to do this in order to be able to work on it.
|
---|
2836 | //
|
---|
2837 | // This version has been put in the repository in order to be
|
---|
2838 | // able to share the further development with others.
|
---|
2839 | //
|
---|
2840 | // If you need something working, wait or take an earlier one.
|
---|
2841 | // See file README.
|
---|
2842 | //
|
---|
2843 | // Revision 1.4 2000/01/25 08:36:23 petry
|
---|
2844 | // The pixelization in previous versions was buggy.
|
---|
2845 | // This is the first version with a correct pixelization.
|
---|
2846 | //
|
---|
2847 | // Revision 1.3 2000/01/20 18:22:17 petry
|
---|
2848 | // Found little bug which makes camera crash if it finds a photon
|
---|
2849 | // of invalid wavelength. This bug is now fixed and the range
|
---|
2850 | // of valid wavelengths extended to 290 - 800 nm.
|
---|
2851 | // This is in preparation for the NSB simulation to come.
|
---|
2852 | // Dirk
|
---|
2853 | //
|
---|
2854 | // Revision 1.2 1999/11/19 08:40:42 harald
|
---|
2855 | // Now it is possible to compile the camera programm under osf1.
|
---|
2856 | //
|
---|
2857 | // Revision 1.1.1.1 1999/11/05 11:59:31 harald
|
---|
2858 | // This the starting point for CVS controlled further developments of the
|
---|
2859 | // camera program. The program was originally written by Jose Carlos.
|
---|
2860 | // But here you can find a "rootified" version to the program. This means
|
---|
2861 | // that there is no hbook stuff in it now. Also the output of the
|
---|
2862 | // program changed to the MagicRawDataFormat.
|
---|
2863 | //
|
---|
2864 | // The "rootification" was done by Dirk Petry and Harald Kornmayer.
|
---|
2865 | //
|
---|
2866 | // Revision 1.3 1999/10/22 15:01:28 petry
|
---|
2867 | // version sent to H.K. and N.M. on Fri Oct 22 1999
|
---|
2868 | //
|
---|
2869 | // Revision 1.2 1999/10/22 09:44:23 petry
|
---|
2870 | // first synthesized version which compiles and runs without crashing;
|
---|
2871 | //
|
---|
2872 | // Revision 1.1.1.1 1999/10/21 16:35:10 petry
|
---|
2873 | // first synthesised version
|
---|
2874 | //
|
---|
2875 | // Revision 1.13 1999/03/15 14:59:05 gonzalez
|
---|
2876 | // camera-1_1
|
---|
2877 | //
|
---|
2878 | // Revision 1.12 1999/03/02 09:56:10 gonzalez
|
---|
2879 | // *** empty log message ***
|
---|
2880 | //
|
---|
2881 | //
|
---|
2882 | //!@}
|
---|
2883 |
|
---|
2884 | //=EOF
|
---|