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

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