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

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