/* ---------------------------------------------------------------------- * * 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) * * ---------------------------------------------------------------------- */ #include #include #include #include #include "version.h" #include "diag.h" #include "init.h" #include "header.h" 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 CerHeader *cheadp; /* var inited in header.c */ extern RflHeader *rheadp; /* var inited in header.c */ /* Prototypes */ void clean(void); FILE *GetNextFile(char *cername); static int GetEvent(void); static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile); void main(void) { long event = 0L; /* event counter */ /* Read init & geometry parms, init files and vars */ init(NULL); /* 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); /* Close file */ Log(RFLF_CLOSE_LOG); fclose(rflfile); /* Clean memory and exit */ clean(); Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION)); } /* end of main */ static int ProcessEvent(CerHeader *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) ); /* Loop on data blocks */ while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile)) && strncmp((char *)Photons, "EVTE", 4)) /* Loop on phs inside block: exit when wlength is 0 */ { for (ph=0; phCorePos[0][0]; Photons[ph].y -= cheadp->CorePos[1][0]; /* Increment number of photons */ phs++; /* Calculate some quantities */ theta = (float) acos(sqrt( 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 */ { /* 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 RflHeader 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(RflHeader), 1, rflfile); /* If there was an overflow, append data from tempfile */ if (overflow) { /* Unload data from CPhotons */ fwrite(CPhotons, sizeof(cphoton), cphs, tmpf); /* Transfer data */ fseek(tmpf, 0L, SEEK_SET); /* Start from the beginning */ while (overflow--) { fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf); fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile); } /* Reload data in CPhotons */ fread (CPhotons, sizeof(cphoton), cphs, tmpf); /* Close (and delete) temp file */ fclose(tmpf); } /* Write (remaining) cphotons */ fwrite(CPhotons, sizeof(cphoton), cphs, rflfile); /* Write "end of event" flag */ fwrite(FLAG_END_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile); Debug(">>> %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 */ 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(CerHeader), 1, cerfile) || strncmp(cheadp->EVTH, "RUNH", 4))) { isWrong = TRUE; } else do { /* Push fileptr in case it is a "EVTH" */ readp = ftell(cerfile); /* Error: exit loop */ if (1 != fread(cheadp, sizeof(CerHeader), 1, cerfile)) { isWrong = TRUE; break; } /* Ok: set found at TRUE and exit loop */ if (strncmp(cheadp->EVTH, "EVTH", 4) == 0 && first_Event <= cheadp->EvtNumber && 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 *f = 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 ((f=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 (f == NULL); /* Loop until a readable file is found */ /* If a new cerfile is met, write "START_OF_RUN" */ if (f) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); return f; } /* 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 */