source: trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c@ 1602

Last change on this file since 1602 was 1535, checked in by bigongia, 22 years ago
Version 0.6. Changed output format: added run header, changed event header, added ascii parameter files attached at the end of every output file to keep all info.
File size: 16.1 KB
Line 
1/* ----------------------------------------------------------------------
2 *
3 * Created: Thu May 7 16:24:22 1998
4 * Author: Jose Carlos Gonzalez
5 * Purpose: Program for reflector simulation
6 * Notes: See files README for details
7 *
8 * ----------------------------------------------------------------------
9 *
10 * Modified: Tue Jan 16 10:10:10 2001
11 * Authors: Denis Bastieri + Ciro Bigongiari
12 * Purpose: Program for reflector simulation (Single file input)
13 *
14 * ----------------------------------------------------------------------
15 *
16 * Modified: October 2002
17 * Author: Abelardo Moralejo
18 * Version 0.6: introduced a run header and changed event header to keep
19 * all information from CORSIKA. In addition, now the magic.def, axisdev.dat
20 * and reflectivity.dat are attached at the end of the reflector output file
21 * to keep also all the parameters with which the program is run.
22 *
23 * ----------------------------------------------------------------------
24 */
25
26#include <stdio.h>
27#include <stdlib.h>
28#include <string.h>
29#include <math.h>
30#include "version.h"
31#include "diag.h"
32#include "init.h"
33#include "header.h"
34
35#define MAX(x,y) ((x)>(y)? (x) : (y))
36#define MIN(x,y) ((x)>(y)? (y) : (x))
37
38FILE *rflfile = NULL, /* reflector (output) file */
39*cerfile = NULL; /* current cerfile */
40
41cphoton *CPhotons = NULL; /* cphoton array */
42static photon Photons[PH_IN_DATABLOCK]; /* photon array */
43
44static long readp = 0L; /* read ptr (on cerfile) */
45static long writep = 0L; /* write ptr (on rflfile) */
46
47extern CerEventHeader *cheadp; /* var inited in header.c */
48extern RflEventHeader *rheadp; /* var inited in header.c */
49
50/* Prototypes */
51void clean(void);
52FILE *GetNextFile(char *cername);
53static int GetEvent(void);
54static int ProcessEvent(CerEventHeader *cheadp, FILE *cerfile, FILE *rflfile);
55void attach(FILE *f1, FILE *f2);
56
57
58FILE *chkf = NULL;
59long myloop;
60
61Event_end* evt_end;
62
63main(void)
64{
65 extern char axisdev_filename[256], reflectivity_filename[256];
66 extern char geo_filename[256];
67 FILE *dummy;
68 long event = 0L; /* event counter */
69
70 /* Read init & geometry parms, init files and vars */
71 init(NULL);
72
73 /* chkf = fopen("check", "w"); */
74
75 /* Processing loop */
76 while(event < max_Events && GetEvent())
77 {
78 if (ProcessEvent(cheadp, cerfile, rflfile))
79 event++;
80 if (event % SHOW_ME == 0)
81 Log(INFO_EVENT_LOG,
82 cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal);
83 }
84
85 /* Writing final flags */
86 fwrite(FLAG_END_OF_FILE, SIZE_OF_FLAGS, 1, rflfile);
87
88 /* Attach parameter files used for running the program: */
89 dummy = fopen (geo_filename,"r");
90 attach (rflfile, dummy);
91 fclose(dummy);
92
93 dummy = fopen (axisdev_filename,"r");
94 attach (rflfile, dummy);
95 fclose(dummy);
96
97 dummy = fopen (reflectivity_filename,"r");
98 attach (rflfile, dummy);
99 fclose(dummy);
100
101 /* Close reflector output file */
102 Log(RFLF_CLOSE_LOG);
103 fclose(rflfile);
104 /* fclose(chkf); */
105
106 /* Clean memory and exit */
107 clean();
108 Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION));
109
110 } /* end of main */
111
112static int ProcessEvent(CerEventHeader *cheadp, FILE *cerfile, FILE *rflfile)
113{ extern int absorption(float wlen, float height, float theta);
114 extern int ph2cph(photon *ph, cphoton *cph);
115 extern float Omega[3][3];
116 extern float OmegaI[3][3];
117 extern float OmegaCT[3][3];
118 extern float OmegaICT[3][3];
119 extern void makeOmega(float theta, float phi);
120 extern void makeOmegaI(float theta, float phi);
121
122
123 /* Various counters: phs = absphs + refphs[0..3] + cphs */
124 long phs, /* Number of incoming photons */
125 absphs, /* Photons absorbed */
126 refphs[4], /* Photons not reflected */
127 cphs; /* Number of cphotons */
128
129 FILE *tmpf=NULL; /* Temp fp to handle o/f */
130 size_t read; /* items read: != 1 on error */
131 int overflow=0, /* how many o/f on cphs */
132 ref_type, /* ret value from reflection */
133 ph; /* photon number */
134
135 float first = 1e8f, /* Photon arrival times */
136 last = 0;
137
138 /* photon quantities */
139 float wlen, /* photon wavelength */
140 theta; /* photon zenith angle */
141
142 /* Reset counters and set fileptr */
143 phs= absphs= refphs[0]= refphs[1]= refphs[2]= refphs[3]= cphs= 0;
144 writep = ftell(rflfile);
145
146 /* Translate the EVTH from cerfile to rflfile */
147 TranslateHeader(rheadp, cheadp);
148
149 /* Calculate OmegaCT matrices */
150 if (is_Fixed_Target)
151 { makeOmega(fixed_Theta, fixed_Phi);
152 makeOmegaI(fixed_Theta, fixed_Phi); }
153 else
154 { makeOmega(cheadp->Theta, cheadp->Phi);
155 makeOmegaI(cheadp->Theta, cheadp->Phi); }
156 memcpy( OmegaCT, Omega, 9*sizeof(float) );
157 memcpy( OmegaICT, OmegaI, 9*sizeof(float) );
158
159 /* ADDED AM May 2002: now we take into account the telescope
160 * position chosen in the Corsika input card via the CERTEL
161 * option, which must be supplied also in the reflector input
162 * card with the option "telescope_position x y". Otherwise
163 * x = y = 0 is assumed, which is the standard mode of
164 * generation for MAGIC.
165 * Here we change rheadp so that the CorePos written to the
166 * .rfl file is referred to the telescope position. However,
167 * I believe that the original "CorePos" vector points
168 * from the core to the origin of coordinates of Corsika, and
169 * therefore after the subtraction of (Telescope_x,Telescope_y)
170 * the resulting vector CorePos points from the core to the
171 * telescope!
172 */
173
174 rheadp->CorePos[0][0] += Telescope_x;
175 rheadp->CorePos[1][0] += Telescope_y;
176
177 /* Loop on data blocks */
178
179 while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile)))
180 /* Loop on phs inside block: exit when wlength is 0 */
181 {
182 /* If "event end" flag is found, read relevant quantities
183 * from event end subblock (added June 2002, AM):
184 */
185
186 if (strncmp((char *)Photons, "EVTE", 4) == 0 )
187 {
188 evt_end = (Event_end*) Photons;
189 rheadp->longi_Nmax = evt_end->longi_Nmax;
190 rheadp->longi_t0 = evt_end->longi_t0;
191 rheadp->longi_tmax = evt_end->longi_tmax;
192 rheadp->longi_a = evt_end->longi_a;
193 rheadp->longi_b = evt_end->longi_b;
194 rheadp->longi_c = evt_end->longi_c;
195 rheadp->longi_chi2 = evt_end->longi_chi2;
196 break;
197 }
198
199 for (ph=0; ph<PH_IN_DATABLOCK; ph++)
200 {
201 /* Added July 2002, AM: check integrity of photon info:
202 Sometimes we found NaN values inside cerXXXXX file.
203 */
204 if (isnan(Photons[ph].w) || isnan(Photons[ph].x) ||
205 isnan(Photons[ph].y) || isnan(Photons[ph].u) ||
206 isnan(Photons[ph].v) || isnan(Photons[ph].t) ||
207 isnan(Photons[ph].h))
208 {
209 printf("Warning: skipped one photon because its data contained Not-a-Number!\n");
210 continue;
211 }
212
213 if (! (Photons[ph].w > 0.) ) break;
214
215 CPhotons[cphs].w = Photons[ph].w;
216 Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.);
217
218 /* AM Nov 2001: we found that sometimes the value
219 stored in Photons[ph].w is not correct, and can result in a
220 wavelength beyond 600 nm, which makes the program crash later.
221
222 AM July 2002: we now simply skip the photon if the wavelength
223 is not in the expected range, which we now take from the corsika
224 event header (just in case we would change it in the future).
225 */
226
227 if (wlen < cheadp->CWaveLower || wlen > cheadp->CWaveUpper)
228 {
229 printf("Warning: skipped one photon with strange wavelength: %f nm\n", wlen);
230 continue;
231 }
232
233 /* ADDED AM May 2002: now we take into account the telescope
234 * position chosen in the Corsika input card via the CERTEL
235 * option, which must be supplied also in the reflector input
236 * card with the option "telescope_position x y". Otherwise
237 * x = y = 0 is assumed, which is the standard mode of
238 * generation for MAGIC.
239 */
240
241 Photons[ph].x -= cheadp->CorePos[0][0]+Telescope_x;
242 Photons[ph].y -= cheadp->CorePos[1][0]+Telescope_y;
243
244 /* Increment number of photons */
245 phs++;
246
247 /* Calculate some quantities */
248 theta = (float) acos(sqrt(
249 MAX(0., 1.f - Photons[ph].u*Photons[ph].u
250 - Photons[ph].v*Photons[ph].v)));
251
252 /* Check absorption */
253
254
255 if (absorption(wlen, Photons[ph].h, theta))
256 absphs++;
257
258 /* Check reflection */
259 else if (0 != (ref_type =
260 ph2cph(&Photons[ph], &CPhotons[cphs])))
261 refphs[ref_type-1]++;
262
263 else /* Photon passed */
264 {
265 Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph,
266 Photons[ph].x,Photons[ph].y,
267 Photons[ph].u,Photons[ph].v,theta);
268 Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x,
269 CPhotons[cphs].y);
270
271
272 /* Update first/last arrival times */
273 if (first > CPhotons[cphs].t) first = CPhotons[cphs].t;
274 if ( last < CPhotons[cphs].t) last = CPhotons[cphs].t;
275
276 /* Update cphs */
277 if (++cphs == NR_OF_CPHOTONS) /* Overflow */
278 {
279 /* if it is the first o/f open a tempfile */
280 if (overflow++ == 0)
281 { Log("Spooling... ");
282 tmpf = tmpfile();
283 if (tmpf == NULL) FatalError(TEMP_ERROR_FTL); }
284
285 /* Write now the cphoton array */
286 fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
287
288 /* Reset the cphs counter */
289 cphs = 0;
290
291 } /* if overflow */
292 } /* else (=Photon passed) */
293 } /* end of loop inside datablock */
294 } /* end of loop on datablocks */
295
296 /* Check if there was an error while reading cerfile */
297 if (read != 1) fseek(rflfile, writep, SEEK_SET);
298
299 else /* no error: write the new event */
300
301 { /* Write "start of event" flag */
302 fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);
303
304 /* Update arrival times */
305 rheadp->TimeFirst = first;
306 rheadp->TimeLast = last;
307
308 /* Update RflEventHeader with info on ph/cph nrs and write it */
309 rheadp->CORSIKAPhs = phs;
310 rheadp->AtmAbsPhs = absphs;
311 rheadp->MirrAbsPhs = refphs[0];
312 rheadp->OutOfMirrPhs = refphs[1];
313 rheadp->BlackSpotPhs = refphs[2];
314 rheadp->OutOfChamPhs = refphs[3];
315 rheadp->CPhotons = (long) overflow * NR_OF_CPHOTONS + cphs;
316 fwrite(rheadp, sizeof(RflEventHeader), 1, rflfile);
317
318/* for (myloop=0; myloop<sizeof(RflEventHeader)/4; myloop++)
319 * fprintf(chkf, "%e ", *((float *)rheadp+myloop));
320 * fputc('\n', chkf);
321 */
322
323 /* If there was an overflow, append data from tempfile */
324 if (overflow)
325 {
326 /* Unload data from CPhotons */
327 fwrite(CPhotons, sizeof(cphoton), cphs, tmpf);
328
329 /* Transfer data */
330 fseek(tmpf, 0L, SEEK_SET); /* Start from the beginning */
331 while (overflow--)
332 {
333 fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
334 fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile);
335 }
336
337 /* Reload data in CPhotons */
338 fread (CPhotons, sizeof(cphoton), cphs, tmpf);
339
340 /* Close (and delete) temp file */
341 fclose(tmpf);
342 }
343
344 /* Write (remaining) cphotons */
345
346 fwrite(CPhotons, sizeof(cphoton), cphs, rflfile);
347
348 /* Write "end of event" flag */
349 fwrite(FLAG_END_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);
350
351 Debug(">>> %6ld = %6ld + %6ld + %6ld + %6ld + %6ld + %6f\n",
352 phs, absphs,
353 refphs[0], refphs[1], refphs[2], refphs[3],
354 rheadp->CPhotons); }
355
356 return read == 1;
357 } /* end of ProcessEvent */
358
359static int GetEvent(void)
360{
361 int found = FALSE, /* event found */
362 isWrong = FALSE; /* cerfile is wrong */
363 static int newFile = TRUE; /* if TRUE, check if cerfile is valid */
364
365 CerRunHeader cRunHead;
366
367 do
368 {
369 /* In case the just-opened file is a valid cerfile,
370 starting with a RUNH, loop until a valid event is found:
371 id est with: first_Event <= EvtNumber <= last_Event
372 and low_Ecut <= Etotal <= high_Ecut
373 If the search was successful, "found" is set to TRUE.
374 If there are reading errors, "isWrong" is set to TRUE. */
375
376 if (newFile
377 && (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile)
378 || strncmp(cheadp->EVTH, "RUNH", 4)))
379 { isWrong = TRUE; }
380
381 else do
382 {
383 if (newFile)
384 {
385 /* Write Reflector "run header" (one per cer file!): */
386 memcpy(&cRunHead, cheadp, sizeof(CerEventHeader));
387 fwrite(&cRunHead, sizeof(CerRunHeader), 1, rflfile);
388 }
389
390 /* Push fileptr in case it is a "EVTH" */
391 readp = ftell(cerfile);
392
393 /* Error: exit loop */
394 if (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile))
395 { isWrong = TRUE; break; }
396
397 /* Ok: set found at TRUE and exit loop */
398 if (strncmp(cheadp->EVTH, "EVTH", 4) == 0
399 && first_Event <= (long)cheadp->EvtNumber
400 && (long)cheadp->EvtNumber <= last_Event
401 && low_Ecut <= cheadp->Etotal
402 && cheadp->Etotal <= high_Ecut)
403 { found = TRUE; }
404
405 /* File is finished: exit loop */
406 else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0)
407 break;
408
409 } while(found == FALSE);
410
411 /* If there was an error, log it */
412 if (isWrong) Log(CERF_ERROR_LOG, cername);
413
414 /* If a new event was not found, get, if any, a new cerfile.
415 "newFile" is set in this case, to check the validity of
416 just-opened cerfile. If GetNextFile return a NULL (no
417 new cerfile), then exit thoroughly. */
418
419 if (found) newFile = FALSE;
420 else
421 { if ((cerfile=GetNextFile(cername)) == NULL) break;
422 newFile = TRUE; }
423
424 } while(found == FALSE);
425
426 return found;
427} /* end of GetEvent */
428
429/* Files are all CER-like files (cerfiles).
430 Filenames are written one per line and may be followed
431 by a ",FIRST_EVT,LAST_EVT" directive, such as:
432 -> cer19999 15 167
433 telling to deal with events from 15 to 167 in file "cer19999"
434 The file list may be appended to the input file or can be stored
435 in a separate file pointed by an '@' directive contained in the
436 input file. In both cases, file referral follows the tag
437 "cer_files". Files not accessible are thoroughly skipped.
438 GetNextFile gets and opens the next readable cerfile. */
439
440FILE *GetNextFile(char *cername)
441{ FILE *inputfile = NULL; /* return value (cerfile ptr) */
442 char *value_ptr; /* ptr at parm value */
443 extern char line[]; /* white chars (init) */
444 extern char whites[]; /* parsing buf. (init) */
445
446 /* Close, if any, the previous cerfile, writing "END_OF_RUN". */
447 if (cerfile)
448 { fclose(cerfile);
449 fwrite(FLAG_END_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); }
450
451 /* Parse next line in filelist */
452 do
453 { if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) break;
454 if ((value_ptr=strtok(line, whites)) == NULL) continue;
455
456 /* If you found a line with some meaning, try to open the file */
457 if ((inputfile=fopen(value_ptr, "r")) == NULL)
458 Message(CERF_NFND__MSG, value_ptr);
459
460 else /* Cerfile is readable */
461 { /* Store its name in "cername" */
462 Log(CERF_OPEN__LOG, strcpy(cername, value_ptr));
463
464 /* Read, if any, event bounds */
465 if ((value_ptr=strtok(NULL, whites)) == NULL)
466 { first_Event = 0;
467 last_Event = 1000000; }
468 else
469 { first_Event = atol(value_ptr);
470 value_ptr = strtok(NULL, whites);
471 last_Event = value_ptr ? atol(value_ptr) : 1000000; }
472
473 /* Check boundaries */
474 if (first_Event > last_Event)
475 { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername);
476 first_Event = 0;
477 last_Event = 1000000; }}}
478
479 while (inputfile == NULL); /* Loop until a readable file is found */
480
481 /* If a new cerfile is met, write "START_OF_RUN" */
482 if (inputfile) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);
483
484 return inputfile;
485} /* end of GetNextFile */
486
487/* This function frees up allocated memory before exiting */
488
489void clean(void)
490{ Message(MEM__FREE__MSG);
491
492 if (ct_data)
493 { free(ct_data); Log(VECT_FREE__LOG, "ct_data"); }
494
495 if (Reflectivity[0])
496 { free(Reflectivity[0]); Log(VECT_FREE__LOG, "Reflectivity[0]"); }
497 if (Reflectivity[1])
498 { free(Reflectivity[1]); Log(VECT_FREE__LOG, "Reflectivity[1]"); }
499
500 if (AxisDeviation[0])
501 { free(AxisDeviation[0]); Log(VECT_FREE__LOG, "AxisDeviation[0]"); }
502 if (AxisDeviation[1])
503 { free(AxisDeviation[1]); Log(VECT_FREE__LOG, "AxisDeviation[1]"); }
504
505 if (ct_Focal)
506 { free(ct_Focal); Log(VECT_FREE__LOG, "ct_Focal"); }
507
508 if (CPhotons)
509 { free(CPhotons); Log(VECT_FREE__LOG, "CPhotons"); }
510} /* end of clean */
511
Note: See TracBrowser for help on using the repository browser.