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

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