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

Last change on this file since 1627 was 1622, checked in by bigongia, 22 years ago
*** empty log message ***
File size: 19.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 */
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 AM Nov 2002: Changed the w range to the fixed values 290 and
292 600 nm. The StarFieldAdder cer files cannot be read otherwise,
293 because the w range is not written in their headers.
294 */
295
296 if (wlen < 290 || wlen > 600)
297 continue;
298
299
300 /* ADDED AM May 2002: now we take into account the telescope
301 * position chosen in the Corsika input card via the CERTEL
302 * option, which must be supplied also in the reflector input
303 * card with the option "telescope_position x y". Otherwise
304 * x = y = 0 is assumed, which is the standard mode of
305 * generation for MAGIC.
306 */
307
308 Photons[ph].x -= cheadp->CorePos[0][0]+Telescope_x;
309 Photons[ph].y -= cheadp->CorePos[1][0]+Telescope_y;
310
311 /* Increment number of photons */
312 phs++;
313
314 /* Calculate some quantities */
315 theta = (float) acos(sqrt(
316 MAX(0., 1.f - Photons[ph].u*Photons[ph].u
317 - Photons[ph].v*Photons[ph].v)));
318
319 /* Check absorption */
320
321
322 if (absorption(wlen, Photons[ph].h, theta))
323 absphs++;
324
325 /* Check reflection */
326 else if (0 != (ref_type =
327 ph2cph(&Photons[ph], &CPhotons[cphs])))
328 refphs[ref_type-1]++;
329
330 else /* Photon passed */
331 {
332 Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph,
333 Photons[ph].x,Photons[ph].y,
334 Photons[ph].u,Photons[ph].v,theta);
335 Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x,
336 CPhotons[cphs].y);
337
338 /* AM Nov 2002: We now count for every event how many
339 * photons were produced by each type.
340 */
341
342 switch (parent_id)
343 {
344 case 2:
345 el_cphs ++;
346 break;
347 case 3:
348 el_cphs ++;
349 break;
350 case 5:
351 mu_cphs ++;
352 break;
353 case 6:
354 mu_cphs ++;
355 break;
356 default:
357 other_cphs++;
358 }
359
360 /* Update first/last arrival times */
361 if (first > CPhotons[cphs].t) first = CPhotons[cphs].t;
362 if ( last < CPhotons[cphs].t) last = CPhotons[cphs].t;
363
364 /* Update cphs */
365 if (++cphs == NR_OF_CPHOTONS) /* Overflow */
366 {
367 /* if it is the first o/f open a tempfile */
368 if (overflow++ == 0)
369 { Log("Spooling... ");
370 tmpf = tmpfile();
371 if (tmpf == NULL) FatalError(TEMP_ERROR_FTL); }
372
373 /* Write now the cphoton array */
374 fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
375
376 /* Reset the cphs counter */
377 cphs = 0;
378
379 } /* if overflow */
380 } /* else (=Photon passed) */
381 } /* end of loop inside datablock */
382 } /* end of loop on datablocks */
383
384 /* Check if there was an error while reading cerfile */
385 if (read != 1) fseek(rflfile, writep, SEEK_SET);
386
387 else /* no error: write the new event */
388
389 { /* Write "start of event" flag */
390 fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);
391
392 /* Update arrival times */
393 rheadp->TimeFirst = first;
394 rheadp->TimeLast = last;
395
396 /* Update RflEventHeader with info on ph/cph nrs and write it */
397 rheadp->CORSIKAPhs = phs;
398 rheadp->AtmAbsPhs = absphs;
399 rheadp->MirrAbsPhs = refphs[0];
400 rheadp->OutOfMirrPhs = refphs[1];
401 rheadp->BlackSpotPhs = refphs[2];
402 rheadp->OutOfChamPhs = refphs[3];
403 rheadp->CPhotons = (long) overflow * NR_OF_CPHOTONS + cphs;
404
405 if (rheadp->CPhotons > 0.)
406 {
407 /* Fill in relative fraction of each particle type producing
408 * Cherenkov light in this even:
409 */
410 rheadp->elec_cph_fraction = (float)el_cphs/(float)rheadp->CPhotons;
411 rheadp->muon_cph_fraction = (float)mu_cphs/(float)rheadp->CPhotons;
412 rheadp->other_cph_fraction = (float)other_cphs/(float)rheadp->CPhotons;
413 }
414 else
415 {
416 rheadp->elec_cph_fraction = 0.;
417 rheadp->muon_cph_fraction = 0.;
418 rheadp->other_cph_fraction = 0.;
419 }
420
421 fwrite(rheadp, sizeof(RflEventHeader), 1, rflfile);
422
423/* for (myloop=0; myloop<sizeof(RflEventHeader)/4; myloop++)
424 * fprintf(chkf, "%e ", *((float *)rheadp+myloop));
425 * fputc('\n', chkf);
426 */
427
428 /* If there was an overflow, append data from tempfile */
429 if (overflow)
430 {
431 /* Unload data from CPhotons */
432 fwrite(CPhotons, sizeof(cphoton), cphs, tmpf);
433
434 /* Transfer data */
435 fseek(tmpf, 0L, SEEK_SET); /* Start from the beginning */
436 while (overflow--)
437 {
438 fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
439 fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile);
440 }
441
442 /* Reload data in CPhotons */
443 fread (CPhotons, sizeof(cphoton), cphs, tmpf);
444
445 /* Close (and delete) temp file */
446 fclose(tmpf);
447 }
448
449 /* Write (remaining) cphotons */
450
451 fwrite(CPhotons, sizeof(cphoton), cphs, rflfile);
452
453 /* Write "end of event" flag */
454 fwrite(FLAG_END_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);
455
456 Debug(">>> %6ld = %6ld + %6ld + %6ld + %6ld + %6ld + %6f\n",
457 phs, absphs,
458 refphs[0], refphs[1], refphs[2], refphs[3],
459 rheadp->CPhotons); }
460
461 return read == 1;
462 } /* end of ProcessEvent */
463
464static int GetEvent(void)
465{
466 int found = FALSE, /* event found */
467 isWrong = FALSE; /* cerfile is wrong */
468 static int newFile = TRUE; /* if TRUE, check if cerfile is valid */
469
470 RflRunHeader RflRunHead;
471
472 extern int atmModel; /* current atm. model */
473
474 do
475 {
476 /* In case the just-opened file is a valid cerfile,
477 starting with a RUNH, loop until a valid event is found:
478 id est with: first_Event <= EvtNumber <= last_Event
479 and low_Ecut <= Etotal <= high_Ecut
480 If the search was successful, "found" is set to TRUE.
481 If there are reading errors, "isWrong" is set to TRUE. */
482
483 if (newFile
484 && (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile)
485 || strncmp(cheadp->EVTH, "RUNH", 4)))
486 { isWrong = TRUE; }
487
488 else do
489 {
490 if (newFile)
491 {
492 /* Write Reflector "run header" (one per cer file!): */
493 memcpy(&RflRunHead, cheadp, sizeof(RflRunHead));
494 RflRunHead.wobble_mode = wobble_position;
495 RflRunHead.atmospheric_model = atmModel;
496 fwrite(&RflRunHead, sizeof(RflRunHead), 1, rflfile);
497 }
498
499 /* Push fileptr in case it is a "EVTH" */
500 readp = ftell(cerfile);
501
502 /* Error: exit loop */
503 if (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile))
504 { isWrong = TRUE; break; }
505
506 /* Ok: set found at TRUE and exit loop */
507 if (strncmp(cheadp->EVTH, "EVTH", 4) == 0
508 && first_Event <= (long)cheadp->EvtNumber
509 && (long)cheadp->EvtNumber <= last_Event
510 && low_Ecut <= cheadp->Etotal
511 && cheadp->Etotal <= high_Ecut)
512 { found = TRUE; }
513
514 /* File is finished: exit loop */
515 else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0)
516 break;
517
518 } while(found == FALSE);
519
520 /* If there was an error, log it */
521 if (isWrong) Log(CERF_ERROR_LOG, cername);
522
523 /* If a new event was not found, get, if any, a new cerfile.
524 "newFile" is set in this case, to check the validity of
525 just-opened cerfile. If GetNextFile return a NULL (no
526 new cerfile), then exit thoroughly. */
527
528 if (found) newFile = FALSE;
529 else
530 { if ((cerfile=GetNextFile(cername)) == NULL) break;
531 newFile = TRUE; }
532
533 } while(found == FALSE);
534
535 return found;
536} /* end of GetEvent */
537
538/* Files are all CER-like files (cerfiles).
539 Filenames are written one per line and may be followed
540 by a ",FIRST_EVT,LAST_EVT" directive, such as:
541 -> cer19999 15 167
542 telling to deal with events from 15 to 167 in file "cer19999"
543 The file list may be appended to the input file or can be stored
544 in a separate file pointed by an '@' directive contained in the
545 input file. In both cases, file referral follows the tag
546 "cer_files". Files not accessible are thoroughly skipped.
547 GetNextFile gets and opens the next readable cerfile. */
548
549FILE *GetNextFile(char *cername)
550{ FILE *inputfile = NULL; /* return value (cerfile ptr) */
551 char *value_ptr; /* ptr at parm value */
552 extern char line[]; /* white chars (init) */
553 extern char whites[]; /* parsing buf. (init) */
554
555 /* Close, if any, the previous cerfile, writing "END_OF_RUN". */
556 if (cerfile)
557 { fclose(cerfile);
558 fwrite(FLAG_END_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); }
559
560 /* Parse next line in filelist */
561 do
562 { if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) break;
563 if ((value_ptr=strtok(line, whites)) == NULL) continue;
564
565 /* If you found a line with some meaning, try to open the file */
566 if ((inputfile=fopen(value_ptr, "r")) == NULL)
567 Message(CERF_NFND__MSG, value_ptr);
568
569 else /* Cerfile is readable */
570 { /* Store its name in "cername" */
571 Log(CERF_OPEN__LOG, strcpy(cername, value_ptr));
572
573 /* Read, if any, event bounds */
574 if ((value_ptr=strtok(NULL, whites)) == NULL)
575 { first_Event = 0;
576 last_Event = 1000000; }
577 else
578 { first_Event = atol(value_ptr);
579 value_ptr = strtok(NULL, whites);
580 last_Event = value_ptr ? atol(value_ptr) : 1000000; }
581
582 /* Check boundaries */
583 if (first_Event > last_Event)
584 { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername);
585 first_Event = 0;
586 last_Event = 1000000; }}}
587
588 while (inputfile == NULL); /* Loop until a readable file is found */
589
590 /* If a new cerfile is met, write "START_OF_RUN" */
591 if (inputfile) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);
592
593 return inputfile;
594} /* end of GetNextFile */
595
596/* This function frees up allocated memory before exiting */
597
598void clean(void)
599{ Message(MEM__FREE__MSG);
600
601 if (ct_data)
602 { free(ct_data); Log(VECT_FREE__LOG, "ct_data"); }
603
604 if (Reflectivity[0])
605 { free(Reflectivity[0]); Log(VECT_FREE__LOG, "Reflectivity[0]"); }
606 if (Reflectivity[1])
607 { free(Reflectivity[1]); Log(VECT_FREE__LOG, "Reflectivity[1]"); }
608
609 if (AxisDeviation[0])
610 { free(AxisDeviation[0]); Log(VECT_FREE__LOG, "AxisDeviation[0]"); }
611 if (AxisDeviation[1])
612 { free(AxisDeviation[1]); Log(VECT_FREE__LOG, "AxisDeviation[1]"); }
613
614 if (CPhotons)
615 { free(CPhotons); Log(VECT_FREE__LOG, "CPhotons"); }
616} /* end of clean */
617
Note: See TracBrowser for help on using the repository browser.