/* ---------------------------------------------------------------------- * * 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 */ extern void SetAtmModel(int model, float ol); /* from atm.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(GetEvent()) { if (event >= max_Events) { fwrite(FLAG_END_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); break; } 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); extern int wobble_position; /* 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 */ /* AM, Nov 2002 */ /* Counters of the number of cphs produced by different particle types: */ long mu_cphs, /* by muons */ el_cphs, /* by electrons */ other_cphs; /* by other particles */ int parent_id; /* Id of the particle producing the C-photon */ 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 */ float tel_theta, tel_phi; /* Indicate telescope orientation, following * Corsika's criterium: theta is the zenith * angle, but phi is the azimuth of the vector * along the telescope axis going from the * camera to the dish, measured from north, * anticlockwise (0-2pi). */ /* 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); if (is_Fixed_Target) { tel_theta = fixed_Theta; tel_phi = fixed_Phi; } else /* Use the direction of the primary if no telescope * orientation is given: */ { tel_theta = cheadp->Theta; tel_phi = cheadp->Phi; } /* AM, Nov 2002: added wobble mode option. It refers * to wobble position +-1, 3 defined in TDAS 01-05 */ switch (wobble_position) { case 1: tel_theta = acos(cos(tel_theta)/1.000024369); tel_phi -= atan(0.006981317/sin(tel_theta)); break; case -1: tel_theta = acos(cos(tel_theta)/1.000024369); tel_phi += atan(0.006981317/sin(tel_theta)); break; case 3: tel_theta += 0.006981317; /* 0.4 degrees */ break; default: break; } /* AM, Nov 2004: in case of wobble mode, change also the telescope * theta and phi in the container to be written to the output file * (rheadp) ! */ rheadp->telescopeTheta = tel_theta; rheadp->telescopePhi = tel_phi; /* Calculate OmegaCT matrices */ /* AM 11/2002 : added PI to 2nd argument in makeOmega, to take * into account that theta and phi from Corsika indicate the * direction of the momentum (hence, downgoing) of the primary * particle, phi measured from positive x axis, anticlockwise, and * theta measured from negative z axis (see Corsika manual). In the * call to makeOmega, the second argument must be the azimuth of the * telescope up-going pointing direction. */ makeOmega(tel_theta, tel_phi+M_PI); makeOmegaI(tel_theta, tel_phi+M_PI); 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; /* Initialize cphoton counters: */ el_cphs = mu_cphs = other_cphs = 0; /* 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; /* AM Nov 2002: now we read the type of particle which produced * the Cherenkov photon : */ parent_id = (int)Photons[ph].w/100000; 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). AM Nov 2002: Changed the w range to the fixed values 290 and 600 nm. The StarFieldAdder cer files cannot be read otherwise, because the w range is not written in their headers. */ if (wlen < 290 || wlen > 600) 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); /* AM Nov 2002: We now count for every event how many * photons were produced by each type. */ switch (parent_id) { case 2: el_cphs ++; break; case 3: el_cphs ++; break; case 5: mu_cphs ++; break; case 6: mu_cphs ++; break; default: other_cphs++; } /* 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; if (rheadp->CPhotons > 0.) { /* Fill in relative fraction of each particle type producing * Cherenkov light in this even: */ rheadp->elec_cph_fraction = (float)el_cphs/(float)rheadp->CPhotons; rheadp->muon_cph_fraction = (float)mu_cphs/(float)rheadp->CPhotons; rheadp->other_cph_fraction = (float)other_cphs/(float)rheadp->CPhotons; } else { rheadp->elec_cph_fraction = 0.; rheadp->muon_cph_fraction = 0.; rheadp->other_cph_fraction = 0.; } 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 */ RflRunHeader RflRunHead; const char *AtmModelNames[]={"ATM_NOATMOSPHERE","ATM_90PERCENT","ATM_CORSIKA"}; extern char atmosphere[256]; /* current atm. model */ int model=0; 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) { while (strcmp(atmosphere, AtmModelNames[model])) if (++model == sizeof(AtmModelNames)/sizeof(AtmModelNames[0])) { model = 0; Error(" *** Atm model \"%s\" not found.\n", atmosphere); break; } /* Write Reflector "run header" (one per cer file!): */ memcpy(&RflRunHead, cheadp, sizeof(RflRunHead)); RflRunHead.wobble_mode = wobble_position; RflRunHead.atmospheric_model = model; fwrite(&RflRunHead, sizeof(RflRunHead), 1, rflfile); /* Set up atmosphere (do some necessary calculations) once we * read in the observation level from the run header: */ Log("Setting atm model \"%s\" and observation level %.0f cm.\n", AtmModelNames[model], RflRunHead.HeightLev[0]); SetAtmModel(model,RflRunHead.HeightLev[0]); } /* 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 (CPhotons) { free(CPhotons); Log(VECT_FREE__LOG, "CPhotons"); } } /* end of clean */