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

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