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

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