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

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