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

Last change on this file since 1621 was 1614, checked in by bigongia, 22 years ago
*** empty log message ***
File size: 18.9 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 extern int wobble_position;
123
124 /* Various counters: phs = absphs + refphs[0..3] + cphs */
125 long phs, /* Number of incoming photons */
126 absphs, /* Photons absorbed */
127 refphs[4], /* Photons not reflected */
128 cphs; /* Number of cphotons */
129
130 /* AM, Nov 2002 */
131 /* Counters of the number of cphs produced by different particle types: */
132 long
133 mu_cphs, /* by muons */
134 el_cphs, /* by electrons */
135 other_cphs; /* by other particles */
136
137
138 int parent_id; /* Id of the particle producing the C-photon */
139
140 FILE *tmpf=NULL; /* Temp fp to handle o/f */
141 size_t read; /* items read: != 1 on error */
142 int overflow=0, /* how many o/f on cphs */
143 ref_type, /* ret value from reflection */
144 ph; /* photon number */
145
146 float first = 1e8f, /* Photon arrival times */
147 last = 0;
148
149 /* photon quantities */
150 float wlen, /* photon wavelength */
151 theta; /* photon zenith angle */
152
153 float tel_theta, tel_phi; /* Indicate telescope orientation, following
154 * Corsika's criterium: theta is the zenith
155 * angle, but phi is the azimuth of the vector
156 * along the telescope axis going from the
157 * camera to the dish, measured from north,
158 * anticlockwise (0-2pi).
159 */
160
161 /* Reset counters and set fileptr */
162 phs= absphs= refphs[0]= refphs[1]= refphs[2]= refphs[3]= cphs= 0;
163 writep = ftell(rflfile);
164
165 /* Translate the EVTH from cerfile to rflfile */
166 TranslateHeader(rheadp, cheadp);
167
168 if (is_Fixed_Target)
169 {
170 tel_theta = fixed_Theta;
171 tel_phi = fixed_Phi;
172 }
173 else /* Use the direction of the primary if no telescope
174 * orientation is given: */
175 {
176 tel_theta = cheadp->Theta;
177 tel_phi = cheadp->Phi;
178 }
179
180 /* AM, Nov 2002: added wobble mode option. It refers
181 * to wobble position +-1, 3 defined in TDAS 01-05 */
182
183 switch (wobble_position)
184 {
185 case 1:
186 tel_theta = acos(cos(tel_theta)/1.000024369);
187 tel_phi -= atan(0.006981317/sin(tel_theta));
188 break;
189 case -1:
190 tel_theta = acos(cos(tel_theta)/1.000024369);
191 tel_phi += atan(0.006981317/sin(tel_theta));
192 break;
193 case 3:
194 tel_theta += 0.006981317; /* 0.4 degrees */
195 break;
196 default:
197 break;
198 }
199
200 /* Calculate OmegaCT matrices */
201 /* AM 11/2002 : added PI to 2nd argument in makeOmega, to take
202 * into account that theta and phi from Corsika indicate the
203 * direction of the momentum (hence, downgoing) of the primary
204 * particle, phi measured from positive x axis, anticlockwise, and
205 * theta measured from negative z axis (see Corsika manual). In the
206 * call to makeOmega, the second argument must be the azimuth of the
207 * telescope up-going pointing direction.
208 */
209
210 makeOmega(tel_theta, tel_phi+M_PI);
211 makeOmegaI(tel_theta, tel_phi+M_PI);
212
213 memcpy( OmegaCT, Omega, 9*sizeof(float) );
214 memcpy( OmegaICT, OmegaI, 9*sizeof(float) );
215
216 /* ADDED AM May 2002: now we take into account the telescope
217 * position chosen in the Corsika input card via the CERTEL
218 * option, which must be supplied also in the reflector input
219 * card with the option "telescope_position x y". Otherwise
220 * x = y = 0 is assumed, which is the standard mode of
221 * generation for MAGIC.
222 * Here we change rheadp so that the CorePos written to the
223 * .rfl file is referred to the telescope position. However,
224 * I believe that the original "CorePos" vector points
225 * from the core to the origin of coordinates of Corsika, and
226 * therefore after the subtraction of (Telescope_x,Telescope_y)
227 * the resulting vector CorePos points from the core to the
228 * telescope!
229 */
230
231 rheadp->CorePos[0][0] += Telescope_x;
232 rheadp->CorePos[1][0] += Telescope_y;
233
234 /* Initialize cphoton counters: */
235 el_cphs = mu_cphs = other_cphs = 0;
236
237 /* Loop on data blocks */
238
239 while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile)))
240 /* Loop on phs inside block: exit when wlength is 0 */
241 {
242 /* If "event end" flag is found, read relevant quantities
243 * from event end subblock (added June 2002, AM):
244 */
245
246 if (strncmp((char *)Photons, "EVTE", 4) == 0 )
247 {
248 evt_end = (Event_end*) Photons;
249 rheadp->longi_Nmax = evt_end->longi_Nmax;
250 rheadp->longi_t0 = evt_end->longi_t0;
251 rheadp->longi_tmax = evt_end->longi_tmax;
252 rheadp->longi_a = evt_end->longi_a;
253 rheadp->longi_b = evt_end->longi_b;
254 rheadp->longi_c = evt_end->longi_c;
255 rheadp->longi_chi2 = evt_end->longi_chi2;
256 break;
257 }
258
259 for (ph=0; ph<PH_IN_DATABLOCK; ph++)
260 {
261 /* Added July 2002, AM: check integrity of photon info:
262 Sometimes we found NaN values inside cerXXXXX file.
263 */
264 if (isnan(Photons[ph].w) || isnan(Photons[ph].x) ||
265 isnan(Photons[ph].y) || isnan(Photons[ph].u) ||
266 isnan(Photons[ph].v) || isnan(Photons[ph].t) ||
267 isnan(Photons[ph].h))
268 {
269 printf("Warning: skipped one photon because its data contained Not-a-Number!\n");
270 continue;
271 }
272
273 if (! (Photons[ph].w > 0.) ) break;
274
275 CPhotons[cphs].w = Photons[ph].w;
276
277 /* AM Nov 2002: now we read the type of particle which produced
278 * the Cherenkov photon :
279 */
280 parent_id = (int)Photons[ph].w/100000;
281 Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.);
282
283 /* AM Nov 2001: we found that sometimes the value
284 stored in Photons[ph].w is not correct, and can result in a
285 wavelength beyond 600 nm, which makes the program crash later.
286
287 AM July 2002: we now simply skip the photon if the wavelength
288 is not in the expected range, which we now take from the corsika
289 event header (just in case we would change it in the future).
290 */
291
292 if (wlen < cheadp->CWaveLower || wlen > cheadp->CWaveUpper)
293 {
294 printf("Warning: skipped one photon with strange wavelength: %f nm\n", wlen);
295 continue;
296 }
297
298 /* ADDED AM May 2002: now we take into account the telescope
299 * position chosen in the Corsika input card via the CERTEL
300 * option, which must be supplied also in the reflector input
301 * card with the option "telescope_position x y". Otherwise
302 * x = y = 0 is assumed, which is the standard mode of
303 * generation for MAGIC.
304 */
305
306 Photons[ph].x -= cheadp->CorePos[0][0]+Telescope_x;
307 Photons[ph].y -= cheadp->CorePos[1][0]+Telescope_y;
308
309 /* Increment number of photons */
310 phs++;
311
312 /* Calculate some quantities */
313 theta = (float) acos(sqrt(
314 MAX(0., 1.f - Photons[ph].u*Photons[ph].u
315 - Photons[ph].v*Photons[ph].v)));
316
317 /* Check absorption */
318
319
320 if (absorption(wlen, Photons[ph].h, theta))
321 absphs++;
322
323 /* Check reflection */
324 else if (0 != (ref_type =
325 ph2cph(&Photons[ph], &CPhotons[cphs])))
326 refphs[ref_type-1]++;
327
328 else /* Photon passed */
329 {
330 Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph,
331 Photons[ph].x,Photons[ph].y,
332 Photons[ph].u,Photons[ph].v,theta);
333 Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x,
334 CPhotons[cphs].y);
335
336 /* AM Nov 2002: We now count for every event how many
337 * photons were produced by each type.
338 */
339
340 switch (parent_id)
341 {
342 case 2:
343 el_cphs ++;
344 break;
345 case 3:
346 el_cphs ++;
347 break;
348 case 5:
349 mu_cphs ++;
350 break;
351 case 6:
352 mu_cphs ++;
353 break;
354 default:
355 other_cphs++;
356 }
357
358 /* Update first/last arrival times */
359 if (first > CPhotons[cphs].t) first = CPhotons[cphs].t;
360 if ( last < CPhotons[cphs].t) last = CPhotons[cphs].t;
361
362 /* Update cphs */
363 if (++cphs == NR_OF_CPHOTONS) /* Overflow */
364 {
365 /* if it is the first o/f open a tempfile */
366 if (overflow++ == 0)
367 { Log("Spooling... ");
368 tmpf = tmpfile();
369 if (tmpf == NULL) FatalError(TEMP_ERROR_FTL); }
370
371 /* Write now the cphoton array */
372 fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
373
374 /* Reset the cphs counter */
375 cphs = 0;
376
377 } /* if overflow */
378 } /* else (=Photon passed) */
379 } /* end of loop inside datablock */
380 } /* end of loop on datablocks */
381
382 /* Check if there was an error while reading cerfile */
383 if (read != 1) fseek(rflfile, writep, SEEK_SET);
384
385 else /* no error: write the new event */
386
387 { /* Write "start of event" flag */
388 fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);
389
390 /* Update arrival times */
391 rheadp->TimeFirst = first;
392 rheadp->TimeLast = last;
393
394 /* Update RflEventHeader with info on ph/cph nrs and write it */
395 rheadp->CORSIKAPhs = phs;
396 rheadp->AtmAbsPhs = absphs;
397 rheadp->MirrAbsPhs = refphs[0];
398 rheadp->OutOfMirrPhs = refphs[1];
399 rheadp->BlackSpotPhs = refphs[2];
400 rheadp->OutOfChamPhs = refphs[3];
401 rheadp->CPhotons = (long) overflow * NR_OF_CPHOTONS + cphs;
402
403 if (rheadp->CPhotons > 0.)
404 {
405 /* Fill in relative fraction of each particle type producing
406 * Cherenkov light in this even:
407 */
408 rheadp->elec_cph_fraction = (float)el_cphs/(float)rheadp->CPhotons;
409 rheadp->muon_cph_fraction = (float)mu_cphs/(float)rheadp->CPhotons;
410 rheadp->other_cph_fraction = (float)other_cphs/(float)rheadp->CPhotons;
411 }
412 else
413 {
414 rheadp->elec_cph_fraction = 0.;
415 rheadp->muon_cph_fraction = 0.;
416 rheadp->other_cph_fraction = 0.;
417 }
418
419 fwrite(rheadp, sizeof(RflEventHeader), 1, rflfile);
420
421/* for (myloop=0; myloop<sizeof(RflEventHeader)/4; myloop++)
422 * fprintf(chkf, "%e ", *((float *)rheadp+myloop));
423 * fputc('\n', chkf);
424 */
425
426 /* If there was an overflow, append data from tempfile */
427 if (overflow)
428 {
429 /* Unload data from CPhotons */
430 fwrite(CPhotons, sizeof(cphoton), cphs, tmpf);
431
432 /* Transfer data */
433 fseek(tmpf, 0L, SEEK_SET); /* Start from the beginning */
434 while (overflow--)
435 {
436 fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
437 fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile);
438 }
439
440 /* Reload data in CPhotons */
441 fread (CPhotons, sizeof(cphoton), cphs, tmpf);
442
443 /* Close (and delete) temp file */
444 fclose(tmpf);
445 }
446
447 /* Write (remaining) cphotons */
448
449 fwrite(CPhotons, sizeof(cphoton), cphs, rflfile);
450
451 /* Write "end of event" flag */
452 fwrite(FLAG_END_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);
453
454 Debug(">>> %6ld = %6ld + %6ld + %6ld + %6ld + %6ld + %6f\n",
455 phs, absphs,
456 refphs[0], refphs[1], refphs[2], refphs[3],
457 rheadp->CPhotons); }
458
459 return read == 1;
460 } /* end of ProcessEvent */
461
462static int GetEvent(void)
463{
464 int found = FALSE, /* event found */
465 isWrong = FALSE; /* cerfile is wrong */
466 static int newFile = TRUE; /* if TRUE, check if cerfile is valid */
467
468 RflRunHeader RflRunHead;
469
470 extern int atmModel; /* current atm. model */
471
472 do
473 {
474 /* In case the just-opened file is a valid cerfile,
475 starting with a RUNH, loop until a valid event is found:
476 id est with: first_Event <= EvtNumber <= last_Event
477 and low_Ecut <= Etotal <= high_Ecut
478 If the search was successful, "found" is set to TRUE.
479 If there are reading errors, "isWrong" is set to TRUE. */
480
481 if (newFile
482 && (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile)
483 || strncmp(cheadp->EVTH, "RUNH", 4)))
484 { isWrong = TRUE; }
485
486 else do
487 {
488 if (newFile)
489 {
490 /* Write Reflector "run header" (one per cer file!): */
491 memcpy(&RflRunHead, cheadp, sizeof(RflRunHead));
492 RflRunHead.wobble_mode = wobble_position;
493 RflRunHead.atmospheric_model = atmModel;
494 fwrite(&RflRunHead, sizeof(RflRunHead), 1, rflfile);
495 }
496
497 /* Push fileptr in case it is a "EVTH" */
498 readp = ftell(cerfile);
499
500 /* Error: exit loop */
501 if (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile))
502 { isWrong = TRUE; break; }
503
504 /* Ok: set found at TRUE and exit loop */
505 if (strncmp(cheadp->EVTH, "EVTH", 4) == 0
506 && first_Event <= (long)cheadp->EvtNumber
507 && (long)cheadp->EvtNumber <= last_Event
508 && low_Ecut <= cheadp->Etotal
509 && cheadp->Etotal <= high_Ecut)
510 { found = TRUE; }
511
512 /* File is finished: exit loop */
513 else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0)
514 break;
515
516 } while(found == FALSE);
517
518 /* If there was an error, log it */
519 if (isWrong) Log(CERF_ERROR_LOG, cername);
520
521 /* If a new event was not found, get, if any, a new cerfile.
522 "newFile" is set in this case, to check the validity of
523 just-opened cerfile. If GetNextFile return a NULL (no
524 new cerfile), then exit thoroughly. */
525
526 if (found) newFile = FALSE;
527 else
528 { if ((cerfile=GetNextFile(cername)) == NULL) break;
529 newFile = TRUE; }
530
531 } while(found == FALSE);
532
533 return found;
534} /* end of GetEvent */
535
536/* Files are all CER-like files (cerfiles).
537 Filenames are written one per line and may be followed
538 by a ",FIRST_EVT,LAST_EVT" directive, such as:
539 -> cer19999 15 167
540 telling to deal with events from 15 to 167 in file "cer19999"
541 The file list may be appended to the input file or can be stored
542 in a separate file pointed by an '@' directive contained in the
543 input file. In both cases, file referral follows the tag
544 "cer_files". Files not accessible are thoroughly skipped.
545 GetNextFile gets and opens the next readable cerfile. */
546
547FILE *GetNextFile(char *cername)
548{ FILE *inputfile = NULL; /* return value (cerfile ptr) */
549 char *value_ptr; /* ptr at parm value */
550 extern char line[]; /* white chars (init) */
551 extern char whites[]; /* parsing buf. (init) */
552
553 /* Close, if any, the previous cerfile, writing "END_OF_RUN". */
554 if (cerfile)
555 { fclose(cerfile);
556 fwrite(FLAG_END_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); }
557
558 /* Parse next line in filelist */
559 do
560 { if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) break;
561 if ((value_ptr=strtok(line, whites)) == NULL) continue;
562
563 /* If you found a line with some meaning, try to open the file */
564 if ((inputfile=fopen(value_ptr, "r")) == NULL)
565 Message(CERF_NFND__MSG, value_ptr);
566
567 else /* Cerfile is readable */
568 { /* Store its name in "cername" */
569 Log(CERF_OPEN__LOG, strcpy(cername, value_ptr));
570
571 /* Read, if any, event bounds */
572 if ((value_ptr=strtok(NULL, whites)) == NULL)
573 { first_Event = 0;
574 last_Event = 1000000; }
575 else
576 { first_Event = atol(value_ptr);
577 value_ptr = strtok(NULL, whites);
578 last_Event = value_ptr ? atol(value_ptr) : 1000000; }
579
580 /* Check boundaries */
581 if (first_Event > last_Event)
582 { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername);
583 first_Event = 0;
584 last_Event = 1000000; }}}
585
586 while (inputfile == NULL); /* Loop until a readable file is found */
587
588 /* If a new cerfile is met, write "START_OF_RUN" */
589 if (inputfile) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);
590
591 return inputfile;
592} /* end of GetNextFile */
593
594/* This function frees up allocated memory before exiting */
595
596void clean(void)
597{ Message(MEM__FREE__MSG);
598
599 if (ct_data)
600 { free(ct_data); Log(VECT_FREE__LOG, "ct_data"); }
601
602 if (Reflectivity[0])
603 { free(Reflectivity[0]); Log(VECT_FREE__LOG, "Reflectivity[0]"); }
604 if (Reflectivity[1])
605 { free(Reflectivity[1]); Log(VECT_FREE__LOG, "Reflectivity[1]"); }
606
607 if (AxisDeviation[0])
608 { free(AxisDeviation[0]); Log(VECT_FREE__LOG, "AxisDeviation[0]"); }
609 if (AxisDeviation[1])
610 { free(AxisDeviation[1]); Log(VECT_FREE__LOG, "AxisDeviation[1]"); }
611
612 if (CPhotons)
613 { free(CPhotons); Log(VECT_FREE__LOG, "CPhotons"); }
614} /* end of clean */
615
Note: See TracBrowser for help on using the repository browser.