/* ---------------------------------------------------------------------- * * Created: Thu May 7 16:24:22 1998 * Author: Jose Carlos Gonzalez * Purpose: Program for reflector simulation * Notes: See files README for details * * ---------------------------------------------------------------------- * * Modified: Tue Jan 16 10:10:10 2001 * Authors: Denis Bastieri + Ciro Bigongiari * Purpose: Program for reflector simulation (Single file input) * * ---------------------------------------------------------------------- * * Modified: October 2002 * Author: Abelardo Moralejo * Version 0.6: introduced a run header and changed event header to keep * all information from CORSIKA. In addition, now the magic.def, axisdev.dat * and reflectivity.dat are attached at the end of the reflector output file * to keep also all the parameters with which the program is run. * * ---------------------------------------------------------------------- */ #include #include #include #include #include "version.h" #include "diag.h" #include "init.h" #include "header.h" #define MAX(x,y) ((x)>(y)? (x) : (y)) #define MIN(x,y) ((x)>(y)? (y) : (x)) FILE *rflfile = NULL, /* reflector (output) file */ *cerfile = NULL; /* current cerfile */ cphoton *CPhotons = NULL; /* cphoton array */ static photon Photons[PH_IN_DATABLOCK]; /* photon array */ static long readp = 0L; /* read ptr (on cerfile) */ static long writep = 0L; /* write ptr (on rflfile) */ extern CerEventHeader *cheadp; /* var inited in header.c */ extern RflEventHeader *rheadp; /* var inited in header.c */ /* Prototypes */ void clean(void); FILE *GetNextFile(char *cername); static int GetEvent(void); static int ProcessEvent(CerEventHeader *cheadp, FILE *cerfile, FILE *rflfile); void attach(FILE *f1, FILE *f2); FILE *chkf = NULL; long myloop; Event_end* evt_end; main(void) { extern char axisdev_filename[256], reflectivity_filename[256]; extern char geo_filename[256]; FILE *dummy; long event = 0L; /* event counter */ /* Read init & geometry parms, init files and vars */ init(NULL); /* chkf = fopen("check", "w"); */ /* Processing loop */ while(event < max_Events && GetEvent()) { if (ProcessEvent(cheadp, cerfile, rflfile)) event++; if (event % SHOW_ME == 0) Log(INFO_EVENT_LOG, cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal); } /* Writing final flags */ fwrite(FLAG_END_OF_FILE, SIZE_OF_FLAGS, 1, rflfile); /* Attach parameter files used for running the program: */ dummy = fopen (geo_filename,"r"); attach (rflfile, dummy); fclose(dummy); dummy = fopen (axisdev_filename,"r"); attach (rflfile, dummy); fclose(dummy); dummy = fopen (reflectivity_filename,"r"); attach (rflfile, dummy); fclose(dummy); /* Close reflector output file */ Log(RFLF_CLOSE_LOG); fclose(rflfile); /* fclose(chkf); */ /* Clean memory and exit */ clean(); Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION)); } /* end of main */ static int ProcessEvent(CerEventHeader *cheadp, FILE *cerfile, FILE *rflfile) { extern int absorption(float wlen, float height, float theta); extern int ph2cph(photon *ph, cphoton *cph); extern float Omega[3][3]; extern float OmegaI[3][3]; extern float OmegaCT[3][3]; extern float OmegaICT[3][3]; extern void makeOmega(float theta, float phi); extern void makeOmegaI(float theta, float phi); /* Various counters: phs = absphs + refphs[0..3] + cphs */ long phs, /* Number of incoming photons */ absphs, /* Photons absorbed */ refphs[4], /* Photons not reflected */ cphs; /* Number of cphotons */ FILE *tmpf=NULL; /* Temp fp to handle o/f */ size_t read; /* items read: != 1 on error */ int overflow=0, /* how many o/f on cphs */ ref_type, /* ret value from reflection */ ph; /* photon number */ float first = 1e8f, /* Photon arrival times */ last = 0; /* photon quantities */ float wlen, /* photon wavelength */ theta; /* photon zenith angle */ /* Reset counters and set fileptr */ phs= absphs= refphs[0]= refphs[1]= refphs[2]= refphs[3]= cphs= 0; writep = ftell(rflfile); /* Translate the EVTH from cerfile to rflfile */ TranslateHeader(rheadp, cheadp); /* Calculate OmegaCT matrices */ if (is_Fixed_Target) { makeOmega(fixed_Theta, fixed_Phi); makeOmegaI(fixed_Theta, fixed_Phi); } else { makeOmega(cheadp->Theta, cheadp->Phi); makeOmegaI(cheadp->Theta, cheadp->Phi); } memcpy( OmegaCT, Omega, 9*sizeof(float) ); memcpy( OmegaICT, OmegaI, 9*sizeof(float) ); /* ADDED AM May 2002: now we take into account the telescope * position chosen in the Corsika input card via the CERTEL * option, which must be supplied also in the reflector input * card with the option "telescope_position x y". Otherwise * x = y = 0 is assumed, which is the standard mode of * generation for MAGIC. * Here we change rheadp so that the CorePos written to the * .rfl file is referred to the telescope position. However, * I believe that the original "CorePos" vector points * from the core to the origin of coordinates of Corsika, and * therefore after the subtraction of (Telescope_x,Telescope_y) * the resulting vector CorePos points from the core to the * telescope! */ rheadp->CorePos[0][0] += Telescope_x; rheadp->CorePos[1][0] += Telescope_y; /* Loop on data blocks */ while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile))) /* Loop on phs inside block: exit when wlength is 0 */ { /* If "event end" flag is found, read relevant quantities * from event end subblock (added June 2002, AM): */ if (strncmp((char *)Photons, "EVTE", 4) == 0 ) { evt_end = (Event_end*) Photons; rheadp->longi_Nmax = evt_end->longi_Nmax; rheadp->longi_t0 = evt_end->longi_t0; rheadp->longi_tmax = evt_end->longi_tmax; rheadp->longi_a = evt_end->longi_a; rheadp->longi_b = evt_end->longi_b; rheadp->longi_c = evt_end->longi_c; rheadp->longi_chi2 = evt_end->longi_chi2; break; } for (ph=0; ph 0.) ) break; CPhotons[cphs].w = Photons[ph].w; Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.); /* AM Nov 2001: we found that sometimes the value stored in Photons[ph].w is not correct, and can result in a wavelength beyond 600 nm, which makes the program crash later. AM July 2002: we now simply skip the photon if the wavelength is not in the expected range, which we now take from the corsika event header (just in case we would change it in the future). */ if (wlen < cheadp->CWaveLower || wlen > cheadp->CWaveUpper) { printf("Warning: skipped one photon with strange wavelength: %f nm\n", wlen); continue; } /* ADDED AM May 2002: now we take into account the telescope * position chosen in the Corsika input card via the CERTEL * option, which must be supplied also in the reflector input * card with the option "telescope_position x y". Otherwise * x = y = 0 is assumed, which is the standard mode of * generation for MAGIC. */ Photons[ph].x -= cheadp->CorePos[0][0]+Telescope_x; Photons[ph].y -= cheadp->CorePos[1][0]+Telescope_y; /* Increment number of photons */ phs++; /* Calculate some quantities */ theta = (float) acos(sqrt( MAX(0., 1.f - Photons[ph].u*Photons[ph].u - Photons[ph].v*Photons[ph].v))); /* Check absorption */ if (absorption(wlen, Photons[ph].h, theta)) absphs++; /* Check reflection */ else if (0 != (ref_type = ph2cph(&Photons[ph], &CPhotons[cphs]))) refphs[ref_type-1]++; else /* Photon passed */ { Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph, Photons[ph].x,Photons[ph].y, Photons[ph].u,Photons[ph].v,theta); Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x, CPhotons[cphs].y); /* Update first/last arrival times */ if (first > CPhotons[cphs].t) first = CPhotons[cphs].t; if ( last < CPhotons[cphs].t) last = CPhotons[cphs].t; /* Update cphs */ if (++cphs == NR_OF_CPHOTONS) /* Overflow */ { /* if it is the first o/f open a tempfile */ if (overflow++ == 0) { Log("Spooling... "); tmpf = tmpfile(); if (tmpf == NULL) FatalError(TEMP_ERROR_FTL); } /* Write now the cphoton array */ fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf); /* Reset the cphs counter */ cphs = 0; } /* if overflow */ } /* else (=Photon passed) */ } /* end of loop inside datablock */ } /* end of loop on datablocks */ /* Check if there was an error while reading cerfile */ if (read != 1) fseek(rflfile, writep, SEEK_SET); else /* no error: write the new event */ { /* Write "start of event" flag */ fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile); /* Update arrival times */ rheadp->TimeFirst = first; rheadp->TimeLast = last; /* Update RflEventHeader with info on ph/cph nrs and write it */ rheadp->CORSIKAPhs = phs; rheadp->AtmAbsPhs = absphs; rheadp->MirrAbsPhs = refphs[0]; rheadp->OutOfMirrPhs = refphs[1]; rheadp->BlackSpotPhs = refphs[2]; rheadp->OutOfChamPhs = refphs[3]; rheadp->CPhotons = (long) overflow * NR_OF_CPHOTONS + cphs; fwrite(rheadp, sizeof(RflEventHeader), 1, rflfile); /* for (myloop=0; myloop>> %6ld = %6ld + %6ld + %6ld + %6ld + %6ld + %6f\n", phs, absphs, refphs[0], refphs[1], refphs[2], refphs[3], rheadp->CPhotons); } return read == 1; } /* end of ProcessEvent */ static int GetEvent(void) { int found = FALSE, /* event found */ isWrong = FALSE; /* cerfile is wrong */ static int newFile = TRUE; /* if TRUE, check if cerfile is valid */ CerRunHeader cRunHead; do { /* In case the just-opened file is a valid cerfile, starting with a RUNH, loop until a valid event is found: id est with: first_Event <= EvtNumber <= last_Event and low_Ecut <= Etotal <= high_Ecut If the search was successful, "found" is set to TRUE. If there are reading errors, "isWrong" is set to TRUE. */ if (newFile && (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile) || strncmp(cheadp->EVTH, "RUNH", 4))) { isWrong = TRUE; } else do { if (newFile) { /* Write Reflector "run header" (one per cer file!): */ memcpy(&cRunHead, cheadp, sizeof(CerEventHeader)); fwrite(&cRunHead, sizeof(CerRunHeader), 1, rflfile); } /* Push fileptr in case it is a "EVTH" */ readp = ftell(cerfile); /* Error: exit loop */ if (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile)) { isWrong = TRUE; break; } /* Ok: set found at TRUE and exit loop */ if (strncmp(cheadp->EVTH, "EVTH", 4) == 0 && first_Event <= (long)cheadp->EvtNumber && (long)cheadp->EvtNumber <= last_Event && low_Ecut <= cheadp->Etotal && cheadp->Etotal <= high_Ecut) { found = TRUE; } /* File is finished: exit loop */ else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0) break; } while(found == FALSE); /* If there was an error, log it */ if (isWrong) Log(CERF_ERROR_LOG, cername); /* If a new event was not found, get, if any, a new cerfile. "newFile" is set in this case, to check the validity of just-opened cerfile. If GetNextFile return a NULL (no new cerfile), then exit thoroughly. */ if (found) newFile = FALSE; else { if ((cerfile=GetNextFile(cername)) == NULL) break; newFile = TRUE; } } while(found == FALSE); return found; } /* end of GetEvent */ /* Files are all CER-like files (cerfiles). Filenames are written one per line and may be followed by a ",FIRST_EVT,LAST_EVT" directive, such as: -> cer19999 15 167 telling to deal with events from 15 to 167 in file "cer19999" The file list may be appended to the input file or can be stored in a separate file pointed by an '@' directive contained in the input file. In both cases, file referral follows the tag "cer_files". Files not accessible are thoroughly skipped. GetNextFile gets and opens the next readable cerfile. */ FILE *GetNextFile(char *cername) { FILE *inputfile = NULL; /* return value (cerfile ptr) */ char *value_ptr; /* ptr at parm value */ extern char line[]; /* white chars (init) */ extern char whites[]; /* parsing buf. (init) */ /* Close, if any, the previous cerfile, writing "END_OF_RUN". */ if (cerfile) { fclose(cerfile); fwrite(FLAG_END_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); } /* Parse next line in filelist */ do { if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) break; if ((value_ptr=strtok(line, whites)) == NULL) continue; /* If you found a line with some meaning, try to open the file */ if ((inputfile=fopen(value_ptr, "r")) == NULL) Message(CERF_NFND__MSG, value_ptr); else /* Cerfile is readable */ { /* Store its name in "cername" */ Log(CERF_OPEN__LOG, strcpy(cername, value_ptr)); /* Read, if any, event bounds */ if ((value_ptr=strtok(NULL, whites)) == NULL) { first_Event = 0; last_Event = 1000000; } else { first_Event = atol(value_ptr); value_ptr = strtok(NULL, whites); last_Event = value_ptr ? atol(value_ptr) : 1000000; } /* Check boundaries */ if (first_Event > last_Event) { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername); first_Event = 0; last_Event = 1000000; }}} while (inputfile == NULL); /* Loop until a readable file is found */ /* If a new cerfile is met, write "START_OF_RUN" */ if (inputfile) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); return inputfile; } /* end of GetNextFile */ /* This function frees up allocated memory before exiting */ void clean(void) { Message(MEM__FREE__MSG); if (ct_data) { free(ct_data); Log(VECT_FREE__LOG, "ct_data"); } if (Reflectivity[0]) { free(Reflectivity[0]); Log(VECT_FREE__LOG, "Reflectivity[0]"); } if (Reflectivity[1]) { free(Reflectivity[1]); Log(VECT_FREE__LOG, "Reflectivity[1]"); } if (AxisDeviation[0]) { free(AxisDeviation[0]); Log(VECT_FREE__LOG, "AxisDeviation[0]"); } if (AxisDeviation[1]) { free(AxisDeviation[1]); Log(VECT_FREE__LOG, "AxisDeviation[1]"); } if (ct_Focal) { free(ct_Focal); Log(VECT_FREE__LOG, "ct_Focal"); } if (CPhotons) { free(CPhotons); Log(VECT_FREE__LOG, "CPhotons"); } } /* end of clean */