Index: trunk/MagicSoft/Simulation/Detector/Starfield/Makefile =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/Makefile (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/Makefile (revision 341) @@ -0,0 +1,219 @@ +################################################################## +## +## Starfield Generator Makefile +## +## $Id: Makefile,v 1.1 2000-01-21 13:36:36 petry Exp $ +## +################################################################## + +# @T \newpage + +# Linux + +INCLUDE = ../include-GENERAL +INCLUDE_COR = ../include-CORSIKA + +CXX = g++ + +INCLUDES = -I${INCLUDE} \ + -I${INCLUDE_COR} + +WARNINGS = -Wall + +CXXFLAGS = ${WARNINGS} ${INCLUDES} -ansi + +LIBS = -L${RANLIBDIR} -L/usr/local/lib -L/usr/lib -lm + + +#------------------------------------------------------------------------------ + +.SILENT: + +.SUFFIXES: .cxx .o + +SRCS = starfield.cxx\ + star.cxx\ + photon.cxx\ + convertcorsika.cxx\ + rand_un_gen.cxx\ + parameters.cxx\ + ${INCLUDE_COR}/COREventHeader.cxx\ + ${INCLUDE_COR}/CORParticle.cxx\ + ${INCLUDE_COR}/CORStatfile.cxx + +HEADERS = starfield.h\ + star.hxx\ + photon.hxx\ + convertcorsika.h\ + parameters.h\ + COREventHeader.hxx\ + CORParticle.hxx\ + CORStatfile.hxx\ + DllImport.h\ + RConfig.h\ + Rtypes.h\ + ranlib.h + +OBJS = starfield.o\ + star.o\ + photon.o\ + convertcorsika.o\ + rand_un_gen.o\ + parameters.o\ + ${INCLUDE_COR}/COREventHeader.o\ + ${INCLUDE_COR}/CORParticle.o\ + ${INCLUDE_COR}/CORStatfile.o + +########################################################### + +all: starfield + +#If you type 'make depend' this will update the dependencies listed below + +depend: + @makedepend $(SRCS) -fMakefile 2> /dev/null + +starfield: $(OBJS) + @echo "Linking..." $(OBJS) + $(CXX) $(CXXFLAGS) $(OBJS) $(LIBS) -o $@ + @echo "done." + +.cxx.o: + @echo "Compiling " $< + $(CXX) $(CXXFLAGS) -c $< -o $@ + +lclean: + @echo "Cleaning..." + @rm -f *.o core + +clean: + @echo "Cleaning..." + @rm -f $(OBJS) core + +mrproper: clean + @echo "Mr.Proper in action . . ." + @rm -f $(PROGRAM) + +ctags: + @echo "Creating CTAGS file . . ." + @ctags -txw $(SRCS) $(HEADERS) > CTAGS + +etags: + @echo "Creating TAGS file . . ." + @/usr/local/bin/etags -C $(SRCS) $(HEADERS) + +listsrc: + @ls -m $(SRCS) $(HEADERS) | sed 's/,//g' + +redo: clean all + +# @endcode + +# DO NOT DELETE THIS LINE -- make depend depends on it. + +starfield.o: starfield.h /usr/include/stdlib.h /usr/include/features.h +starfield.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h +starfield.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +starfield.o: /usr/include/sys/types.h /usr/include/gnu/types.h +starfield.o: /usr/include/time.h /usr/include/endian.h /usr/include/bytesex.h +starfield.o: /usr/include/sys/select.h /usr/include/selectbits.h +starfield.o: /usr/include/alloca.h /usr/include/string.h /usr/include/math.h +starfield.o: /usr/include/huge_val.h /usr/include/mathcalls.h +starfield.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/float.h +starfield.o: /usr/include/dirent.h /usr/include/direntry.h +starfield.o: /usr/include/posix1_lim.h /usr/include/local_lim.h +starfield.o: /usr/include/linux/limits.h /usr/include/unistd.h +starfield.o: /usr/include/posix_opt.h /usr/include/confname.h +starfield.o: convertcorsika.h photon.hxx star.hxx /usr/include/stdio.h +starfield.o: /usr/include/libio.h /usr/include/_G_config.h +starfield.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stdarg.h +starfield.o: /usr/include/stdio_lim.h parameters.h +star.o: star.hxx /usr/include/math.h /usr/include/features.h +star.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h +star.o: /usr/include/huge_val.h /usr/include/mathcalls.h +star.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/float.h +star.o: /usr/include/string.h +star.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +star.o: /usr/include/stdio.h /usr/include/libio.h /usr/include/_G_config.h +star.o: /usr/include/gnu/types.h +star.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stdarg.h +star.o: /usr/include/stdio_lim.h +photon.o: photon.hxx /usr/include/math.h /usr/include/features.h +photon.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h +photon.o: /usr/include/huge_val.h /usr/include/mathcalls.h +photon.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/float.h +photon.o: /usr/include/stdlib.h +photon.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +photon.o: /usr/include/sys/types.h /usr/include/gnu/types.h +photon.o: /usr/include/time.h /usr/include/endian.h /usr/include/bytesex.h +photon.o: /usr/include/sys/select.h /usr/include/selectbits.h +photon.o: /usr/include/alloca.h +convertcorsika.o: convertcorsika.h /usr/include/string.h +convertcorsika.o: /usr/include/features.h /usr/include/sys/cdefs.h +convertcorsika.o: /usr/include/gnu/stubs.h +convertcorsika.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +convertcorsika.o: photon.hxx /usr/include/math.h /usr/include/huge_val.h +convertcorsika.o: /usr/include/mathcalls.h +convertcorsika.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/float.h +convertcorsika.o: /usr/include/stdlib.h /usr/include/sys/types.h +convertcorsika.o: /usr/include/gnu/types.h /usr/include/time.h +convertcorsika.o: /usr/include/endian.h /usr/include/bytesex.h +convertcorsika.o: /usr/include/sys/select.h /usr/include/selectbits.h +convertcorsika.o: /usr/include/alloca.h +parameters.o: parameters.h /usr/include/stdio.h /usr/include/libio.h +parameters.o: /usr/include/features.h /usr/include/sys/cdefs.h +parameters.o: /usr/include/gnu/stubs.h /usr/include/_G_config.h +parameters.o: /usr/include/gnu/types.h +parameters.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +parameters.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stdarg.h +parameters.o: /usr/include/stdio_lim.h +../include-CORSIKA/COREventHeader.o: ../include-CORSIKA/COREventHeader.hxx +../include-CORSIKA/COREventHeader.o: /usr/include/stdlib.h +../include-CORSIKA/COREventHeader.o: /usr/include/features.h +../include-CORSIKA/COREventHeader.o: /usr/include/sys/cdefs.h +../include-CORSIKA/COREventHeader.o: /usr/include/gnu/stubs.h +../include-CORSIKA/COREventHeader.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +../include-CORSIKA/COREventHeader.o: /usr/include/sys/types.h +../include-CORSIKA/COREventHeader.o: /usr/include/gnu/types.h +../include-CORSIKA/COREventHeader.o: /usr/include/time.h +../include-CORSIKA/COREventHeader.o: /usr/include/endian.h +../include-CORSIKA/COREventHeader.o: /usr/include/bytesex.h +../include-CORSIKA/COREventHeader.o: /usr/include/sys/select.h +../include-CORSIKA/COREventHeader.o: /usr/include/selectbits.h +../include-CORSIKA/COREventHeader.o: /usr/include/alloca.h +../include-CORSIKA/COREventHeader.o: /usr/include/math.h +../include-CORSIKA/COREventHeader.o: /usr/include/huge_val.h +../include-CORSIKA/COREventHeader.o: /usr/include/mathcalls.h +../include-CORSIKA/COREventHeader.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/float.h +../include-CORSIKA/CORParticle.o: ../include-CORSIKA/CORParticle.hxx +../include-CORSIKA/CORParticle.o: /usr/include/stdlib.h +../include-CORSIKA/CORParticle.o: /usr/include/features.h +../include-CORSIKA/CORParticle.o: /usr/include/sys/cdefs.h +../include-CORSIKA/CORParticle.o: /usr/include/gnu/stubs.h +../include-CORSIKA/CORParticle.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +../include-CORSIKA/CORParticle.o: /usr/include/sys/types.h +../include-CORSIKA/CORParticle.o: /usr/include/gnu/types.h +../include-CORSIKA/CORParticle.o: /usr/include/time.h /usr/include/endian.h +../include-CORSIKA/CORParticle.o: /usr/include/bytesex.h +../include-CORSIKA/CORParticle.o: /usr/include/sys/select.h +../include-CORSIKA/CORParticle.o: /usr/include/selectbits.h +../include-CORSIKA/CORParticle.o: /usr/include/alloca.h /usr/include/math.h +../include-CORSIKA/CORParticle.o: /usr/include/huge_val.h +../include-CORSIKA/CORParticle.o: /usr/include/mathcalls.h +../include-CORSIKA/CORParticle.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/float.h +../include-CORSIKA/CORStatfile.o: ../include-CORSIKA/CORStatfile.hxx +../include-CORSIKA/CORStatfile.o: /usr/include/stdlib.h +../include-CORSIKA/CORStatfile.o: /usr/include/features.h +../include-CORSIKA/CORStatfile.o: /usr/include/sys/cdefs.h +../include-CORSIKA/CORStatfile.o: /usr/include/gnu/stubs.h +../include-CORSIKA/CORStatfile.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/stddef.h +../include-CORSIKA/CORStatfile.o: /usr/include/sys/types.h +../include-CORSIKA/CORStatfile.o: /usr/include/gnu/types.h +../include-CORSIKA/CORStatfile.o: /usr/include/time.h /usr/include/endian.h +../include-CORSIKA/CORStatfile.o: /usr/include/bytesex.h +../include-CORSIKA/CORStatfile.o: /usr/include/sys/select.h +../include-CORSIKA/CORStatfile.o: /usr/include/selectbits.h +../include-CORSIKA/CORStatfile.o: /usr/include/alloca.h /usr/include/math.h +../include-CORSIKA/CORStatfile.o: /usr/include/huge_val.h +../include-CORSIKA/CORStatfile.o: /usr/include/mathcalls.h +../include-CORSIKA/CORStatfile.o: /usr/lib/gcc-lib/i386-redhat-linux/2.7.2.3/include/float.h Index: trunk/MagicSoft/Simulation/Detector/Starfield/README =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/README (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/README (revision 341) @@ -0,0 +1,205 @@ +$Id: README,v 1.1 2000-01-21 13:36:36 petry Exp $ + +STARFIELD README + +D. Petry, IFAE, Campus UAB, Spain + + +This is the first useful version of the starfield generator (SG). + +Purpose: + +The SG is needed to simulate the non-diffuse part of the night +sky background (NSB) in images taken by Cherenkov telescopes. +It reads data from star catalogues (presently only SKY2000, see below) +and calculated from the given position of the telescopes optical +axis (in celestial coordinates), the positional and spectral data +for each star from the catalog, how many photons of which wavelength +(presently 290 nm - 800 nm are simulated) will hit a circular surface +of given radius ("mirror radius") in a given time ("integration time"). +For each photon the director cosines are calculated. +Arrival time, wavelength and ground position are randomized within the +integration time, the four wavebands (U, B, V, R) and the mirror area +respectively with flat distributions. + +SG finally writes the generated photons into a binary file of +CORSIKA format which can be read by the "Reflector" simulation. +The output form the Reflector can then be fed to the camera. + +The camera will pixelize the photons and calculate the average +NSB photon rate for each pixel. From this, the camera can then +generate the NSB contribution in each shower image using a Poisson +random generator (at present, 21-1-2000, this is not yet implemented +but it will be very soon) and also taking into account the +atmospheric extinction and its dependence on the zenith angle. + + +Installation and first test: + +1) Adjust the Makefile to your system and type + + make depend + make + + The result should be an executable named "starfield" + +2) Download the star catalog: + + ftp cdsarc.u-strasbg.fr + + (log in as user anonymous) + + cd pub/cats/V/102 + prompt + hash + bin + mget sky*.dat.gz + bye + +3) Unzip the 24 files you have downloaded and put them into some permanent + directory e.g. the Data directory of the Monte Carlo source code. + +4) Edit the parameter file. The distribution includes an + example file named "starfield.par" which looks like this: + + Starfield Generator Parameters, Date: 21-1-2000, Comment: Example (Crab Nebula) + Center of FOV: ira_hours ira_min ira_sec idec_degrees idec_arcmin dec_arcsec: + 05 34 32 +22 00 52.1 + Radius of the FOV for the catalog readout (degrees): + 2.0 + Integration time for the calculation of the number of photons (seconds): + 50e-9 + Mirror radius for the generation of random impact points (meters): + 10.0 + Path inside which the star catalog data can be found: + /usr/users/xf/stardata + Verbosity level (0 = not verbose, 1 = verbose, 2 = very verbose, 3 = very very ...): + 0 + + Note that the there is a header line followed by pairs of an explanatory + line and a data line. You may write in the header line and the explanaroy + lines whatever you like, but the data lines matter. The lines may be as long + as you like but don't add any carriage returns. + + For a first test run, you will only have to edit the path for the + star data which is the one defined by yourself in step (3). + + +5) Test run. + Run the program by entering + + starfield + + The diagnostic output to the screen should be ending with + + Opened starfield.par for reading ... + Starfield Generator Parameters, Date: 21-1-2000, Comment: Example (Crab Nebula) + Position RA DEC: 5 34 32 22 0 52.1 + FOV Radius:2 degrees + Integration Time:5e-08 s + Mirror Radius:10 m + Catalog Data Path: ... your path ... + Verbosity: 0 + SKY2000 - Master Star Catalog - Star Catalog Database, Version 2 + Sande C.B., Warren W.H.Jr., Tracewell D.A., Home A.T., Miller A.C. + + Opened file ...your path.../sky04.dat for reading ... + EOF reached; accepted 0 stars from this segment. + Opened file ...your path.../sky05.dat for reading ... + Warning: star no. 53701440 is bright (Vmag =3, Bmag = 2.85) + and has no Umag measurement. Estimated Umag is 2.2395 + EOF reached; accepted 120 stars from this segment. + Opened file ...your path .../sky06.dat for reading ... + EOF reached; accepted 0 stars from this segment. + Accepted 120 stars in total. + Writing binary Cherenkov file ./cer050220 ... + Done. + Writing binary statistics file ./sta050220 ... + Done. + +6) Using the output. + The two output files (in this case cer050220 and sta050220) are + of the same format as the file for a single event in the CORSIKA + shower simulation. The number in the name is generated from the + RA (h) and DEC (deg) of the telescope position in order to allow + a distinction. + + The further processing has to be done with the reflector. + A sample parameter file for the reflector is the following: + + reflector 0.3 -*- sh -*- + # + # Sample parameters file + # + verbose_level 2 + # + fixed_target 0. 0. + ct_file ../Data/magic.def + output_file starfield.rfl + atm_model ATM_NOATMOSPHERE + data_paths 1 + ... the diretory in which the output of starfield is found + end_file + + Note that the line "fixed_target 0. 0." is a must. + Also the line "atm_model ATM_NOATMOSPHERE" is necessary + because the extinction will be simulated in the camera (see above). + +7) Notes on the input parameters of starfield + + a) Center of FOV: ira_hours ira_min ira_sec idec_degrees idec_arcmin dec_arcsec: + + Here you can put any valid celestial coordinates on the sky. + The program doesn't veto invalid coordinates yet, but the result + is undefined. + + b) Radius of the FOV for the catalog readout (degrees): + + This should be a number larger than the outermost radius of the + CT's field of view (FOV), but not too much larger because there + is a limit to the number of photons which can be stored by the + program and the number of generated photons grows with the number + of stars in the FOV. If you hit the limit, try decreasing the + integration time. + + c) Integration time for the calculation of the number of photons (seconds): + + Since we are dealing with a ray-tracing approach, we need actual + individual photons, not just rates. This integration time is needed + in order to normalize the numbers of generated photons for each star. + The camera program will later calculate the photon rates for each + pixel by dividing by this number. + + Not that this number is given in seconds, but only up to about 100 ns + make sense, i.e. you have to write something like 100e-9 + + d) Mirror radius for the generation of random impact points (meters): + + This number is needed for the flux normalization and the generation + of the ground impact points of the photons. Adjust it to the size + of the simulated telescope. + + e) Path inside which the star catalog data can be found: + + This is the directory containing the SKY2000 data (see point 3) + + f) Verbosity level (0 = not verbose, 1 = verbose, 2 = very verbose, 3 = very very ...): + + Set this value to 0 for normal usage. Increase it to discover + what is going on in detail. + +8) parameter file name + + starfield accepts as its only command line argument the name + of the parameter file. Typing + + starfield my.par + + will make starfield look for the parameter file "my.par" in + the current directory. If no argument is given, it assumes the + name of the parameter file is "starfield.par". + + + + + Index: trunk/MagicSoft/Simulation/Detector/Starfield/convertcorsika.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/convertcorsika.cxx (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/convertcorsika.cxx (revision 341) @@ -0,0 +1,126 @@ +////////////////////////////////////////////////////////////////////////////////////// +// Write photon data in binary corsika-like format. +// Output will be used as an input for the reflector as if it was a shower event. +// +////////////////////////////////////////////////////////////////////////////////////// + +#include "convertcorsika.h" + +COREventHeader cerevth; +CORParticle cerphot; +CORStatfile cerstat; + +int convertcorsika(int id, int photnum, photon phot[], float inttime_s, int verbose){ + + int i,filenum; + float x_cm, y_cm, u, v, lambda_nm, t_ns; + char cor_dir[60]; + char stat_dir[60]; + char cor_file[15]; + char stat_file[15]; + + filenum=id; + + // Event corsika file + + //File labeling. + + strcpy (cor_dir, "./"); + + sprintf(cor_file, "cer%06d",filenum); + + strcat(cor_dir, cor_file); + + + // Fill the header of the corsika-like event. + //The fields in evt.fill are: event number, primary identifier, total energy (GeV), + // first target identifier, first z interaction(cm) , momentum x, momentum y, momentum z, + // zenith angle, azimuth angle, coreposition x (cm), coreposition y(cm) + + ofstream cerfile; + + if(!cerfile) { + cerr<<"Cannot create cerenkov file.\n"; + exit(1); + } + + cerfile.open(cor_dir, ios::out); + + //____________________________________________________________________________________ + + cout<<"Writing binary Cherenkov file"<<" "< + +#include "photon.hxx" +//#include "getphotons.h" + +//include classes for corsika filling events from reflector mc + + +#include "COREventHeader.hxx" +#include "CORParticle.hxx" +#include "CORStatfile.hxx" + + +//File definition + +const char CORSIKA_FILE[15] = "cer999"; +const char STAT_FILE[15]= "sta999"; + + +#ifndef PI +#define PI 3.1415927 +#endif + +//Function definition +int convertcorsika(int id, int photnum, photon phot[], float inttime_s, int verbose); + +#endif + + + + + + + + + + Index: trunk/MagicSoft/Simulation/Detector/Starfield/parameters.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/parameters.cxx (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/parameters.cxx (revision 341) @@ -0,0 +1,124 @@ +#include "parameters.h" + +parameters::parameters(){ // constructor (put invalid values) + ct_ra_h = -999.; + ct_dec_deg = -999.; + ct_ra_rad = -999.; + ct_dec_rad = -999.; + catalog_fov_deg = -999.; + integtime_s = -999.; + mirr_radius_m = -999.; + sprintf(datapath, "."); + verbose = 0; +} + +int parameters::readparameters(ifstream *in){ // read the parameter file + int ira_hours, ira_min, ira_sec; + int idec_degrees, idec_arcmin; + float dec_arcsec; + char dummy[160]; + + in->ipfx(1); // tell the io system that we are reading formatted + + in->getline(dummy, sizeof(dummy), '\n'); + cout << dummy << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + sscanf(dummy, "%d %d %d %d %d %f", &ira_hours , &ira_min, &ira_sec, + &idec_degrees, &idec_arcmin, &dec_arcsec); + cout << "Position RA DEC: "<< ira_hours << " " << ira_min << " " << ira_sec << " " << + idec_degrees << " " << idec_arcmin << " " << dec_arcsec << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + sscanf(dummy, "%f", &catalog_fov_deg); + cout << "FOV Radius:" << catalog_fov_deg << " degrees\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + sscanf(dummy, "%f", &integtime_s); + cout << "Integration Time:" << integtime_s << " s\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + sscanf(dummy, "%f", &mirr_radius_m); + cout << "Mirror Radius:" << mirr_radius_m << " m\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + sscanf( dummy, "%s", datapath); + cout << "Catalog Data Path: " << datapath << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + if(in->eof()) return(in->good()); + + in->getline(dummy, sizeof(dummy), '\n'); + if(verbose) cout << dummy << "\n"; + sscanf( dummy, "%d", &verbose); + cout << "Verbosity: " << verbose << "\n"; + + ct_ra_h = ira_hours + ira_min/60. + ira_sec/3600.; + ct_dec_deg = idec_degrees + idec_arcmin/60. + dec_arcsec/3600.; + ct_ra_rad = ct_ra_h * PI / 12.; + ct_dec_rad = ct_dec_deg * PI / 180.; + + return(in->good()); +} + +void parameters::usage(ostream *out){ + *out << "Starfield Generator Parameters, Date: ......, Comment: ....\n"; + *out << "Center of FOV: ira_hours ira_min ira_sec idec_degrees idec_arcmin dec_arcsec:\n"; + *out << "05 34 32 +22 00 52.1 \n"; + *out << "Radius of the FOV for the catalog readout (degrees):\n"; + *out << "2.0\n"; + *out << "Integration time for the calculation of the number of photons (seconds):\n"; + *out << "50e-9\n"; + *out << "Mirror radius for the generation of random impact points (meters):\n"; + *out << "10.0\n"; + *out << "Path inside which the star catalog data can be found:\n"; + *out << "/usr/users/xy/stardata\n"; + *out << "Verbosity level (0 = not verbose, 1 = verbose, 2 = very verbose, 3 = very very ...):\n"; + *out << "0\n"; + return; +} + +int parameters::printparameters(ostream *out){ // write the parameters to the stream + out->setf(ios::adjustfield); + out->width(PRINTWIDTH); + *out << "optical axis RA (h)" << ct_ra_h; out->width(PRINTWIDTH); + *out << "optical axis DEC (deg)" << ct_dec_deg; out->width(PRINTWIDTH); + *out << "\n"; + return(0); +} +// reminder: +// some of the things from star.h should go into the parameters file + + Index: trunk/MagicSoft/Simulation/Detector/Starfield/parameters.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/parameters.h (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/parameters.h (revision 341) @@ -0,0 +1,40 @@ +#ifndef _PARAMETERS_H_ +#define _PARAMETERS_H_ + +#include +#include +#include + +#ifndef PI +#define PI 3.1415927 +#endif + +#ifndef PRINTWIDTH +#define PRINTWIDTH 12 +#endif + + +class parameters { + + public: + + float ct_ra_h, ct_dec_deg, ct_ra_rad, ct_dec_rad; + float catalog_fov_deg; // Radius of the catalog section to be read in + float integtime_s; // integration time for the calculation of the number of photons + float mirr_radius_m; // mirror radius inside which photons are generated + char datapath[160]; // path inside which the catalog and extinction files are found + int verbose; // verbose flag for switching on/off diagnostic output + + parameters(); // constructor (put invalid values) + int readparameters(ifstream *in); // read the parameter file + void usage(ostream *out); // write a sample parameter file to the stream + int printparameters(ostream *out); // write the parameters to the stream + +}; + +#endif + + + + + Index: trunk/MagicSoft/Simulation/Detector/Starfield/photon.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/photon.cxx (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/photon.cxx (revision 341) @@ -0,0 +1,37 @@ +//////////////////////////////////////////////////////////////////////////////////////////////// +// +// Here we define the functions which will give a random coordinate pair (x,y) in the ground +// for every photon and another function which will give a random arrival time inside a time +// window for each arriving photon. +// +// Mireia Dosil 24-II-99 IFAE +//////////////////////////////////////////////////////////////////////////////////////////////// + +#include "photon.hxx" + + +// Constructor function. + +photon::photon(){ + + u=-999.; + v=-999.; + lambda_nm=-999.; + starnum=-999; + arrtime_sec=0.0; + x_m=0.0; + y_m=0.0; +} + + + + + + + + + + + + + Index: trunk/MagicSoft/Simulation/Detector/Starfield/photon.hxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/photon.hxx (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/photon.hxx (revision 341) @@ -0,0 +1,39 @@ +//////////////////////////////////////////////////////////////////// +// Declaration of the class photon which will be inherited from the +// class star. +///////////////////////////////////////////////////////////////////////// + + +#ifndef _PHOTON_H_ +#define _PHOTON_H_ + +#include +#include +#include + + +class photon { + +public: + + float u; + float v; + int starnum; + float lambda_nm; + float arrtime_sec; + float x_m; + float y_m; + + + photon(); //Constructor function. + +}; + + +#endif + + + + + + Index: trunk/MagicSoft/Simulation/Detector/Starfield/rand_un_gen.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/rand_un_gen.cxx (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/rand_un_gen.cxx (revision 341) @@ -0,0 +1,92 @@ +///////////////////////////////////////////////////////////////////////////////////////////////// +// This function is from 'Numerical Recipes in C: The Art of Scientific Computing +// William H. Press, Brian P. Flannery, Saul A. Teulosky, William T. Vetteling; New York: +// Cambridge University Press, 1990.' +// +// It returns uniform random numbers between 0 and 1. +// +////////////////////////////////////////////////////////////////////////////////////////////////// + +#define MBIG 1000000000 +#define MSEED 161803398 +#define MZ 0 +#define FAC (1.0/MBIG) + + +// According to Knuth, any large MBIG, and any smaller ( but still large ) MSEED can be substituted for the above values. + +float rand_un_gen(long *idum_local) + + //Returns a uniform random deviate between 0.0 and 1.0. Set idum_local to any negative value to initialize or reanitialize the sequence. + +{ + + + static int inext, inextp; + static long ma[56]; //The value 56 is special and should not be modified + static int iff=0; + long mj, mk; + int i, ii, k; + + + if (*idum_local < 0 || iff == 0) { //Initialization + iff=1; + mj=MSEED-(*idum_local > 0 ? -*idum_local : *idum_local); // Initialize ma[55] using the seed idum_local and the large number MSEED. + mj %= MBIG; + ma[55]=mj; + mk=1; + + //Now initialize the rest of the table, in a slightly random order, with numbers that are not specially random. + + + + for (i=1; i<=54; i++) { + ii=(21*i) % 55; + ma[ii]=mk; + mk=mj-mk; + if ( mk < MZ ) mk += MBIG; + mj=ma[ii]; + } + + for (k=1;k<=4;k++) + for(i=1;i<=55;i++){ + ma[i] -= ma[1+(i+30) % 5]; + if(ma[i] 2) fprintf(stdout, "%s\n", catline); + + pos = catline + 27; + strncpy(c8, pos, 8); + sscanf(c8, "%d", &icatnum); + + pos = catline + 118; + strncpy(c2, pos, 2); + sscanf(c2, "%d", &ira_hours); + + pos = catline + 120; + strncpy(c2, pos, 2); + sscanf(c2, "%d", &ira_min); + + pos = catline + 122; + strncpy(c7, pos, 7); + sscanf(c7, "%f", &ra_sec); + + pos = catline + 129; + strncpy(c3, pos, 3); + if( c3[1] == ' ' ){ + c3[1] = '0'; + } + if( c3[2] == ' ' ){ + c3[2] = '0'; + } + sscanf(c3, "%d", &idec_degrees); + + pos = catline + 132; + strncpy(c2, pos, 2); + sscanf(c2, "%d", &idec_arcmin); + + pos = catline + 134; + strncpy(c6, pos, 6); + sscanf(c6, "%f", &dec_arcsec); + + pos = catline + 231; + strncpy(c6, pos, 6); + sscanf(c6, "%f", &vmag); + + pos = catline + 251; + strncpy(c6, pos, 6); + sscanf(c6, "%f", &bmag); + + pos = catline + 271; + strncpy(c6, pos, 6); + sscanf(c6, "%f", &umag); + + ra_h = ira_hours + ira_min/60. + ra_sec/3600.; + dec_deg = idec_degrees + idec_arcmin/60. + dec_arcsec/3600.; + ra_rad = ra_h * PI / 12.; + dec_rad = dec_deg * PI /180.; + + if (verbose > 2) fprintf(stdout, "extracted: %d %d %d %f %d %d %f %f %f %f\n", + icatnum, ira_hours, ira_min, ra_sec, + idec_degrees, idec_arcmin, dec_arcsec, + umag, bmag, vmag); + + return(TRUE); +} + +int star::printstar(){ // write one star's parameters + fprintf(stdout, "%d %f %f %f %f %f %f\n", + icatnum, ra_h, dec_deg, umag, bmag, vmag, rmag); + return(0); +} + +//---------------------------------------------------------------------------- +// @name calcmissingmags +// +// @desc calculate the magnitudes for those wavebands in which no data +// @desc is available assuming a black body and using the V and B mags +// +//---------------------------------------------------------------------------- + +float star::calcmissingmags(int verbose) { // returns effective temperature; -1. = not possible + + float temp; + float tprime; + float nu1_Hz, nu2_Hz, bflux, vflux, rflux, xmag; + + + + if(vmag > -100.){ // valid vmag available + if(bmag < -100.){ // no valid bmag available + cout << "Warning: star no. " << icatnum << + " has no Bmag measurement. Using Bmag = Vmag = "< -0.2){ + temp = 8540. / ( (bmag-vmag) + 0.865 ); + if (verbose > 1) cout << "Star temperature from B-V: T = " << temp << "K\n"; + } + else{ + temp = 12000.; + if (verbose > 1) cout << "Star temperature from B-V: T > " << temp << "K\n"; + } + + // calculate an effective temperature for the Rmag calculation + // tprime = T * k / h + + nu1_Hz = LIGHTSPEED_mps/((VLMIN_nm+VLMAX_nm)/2.*1e-9); + nu2_Hz = LIGHTSPEED_mps/((BLMAX_nm+BLMAX_nm)/2.*1e-9); + vflux = pow(10.,-0.4*vmag-22.42); + bflux = pow(10.,-0.4*bmag-22.42); + tprime = (nu2_Hz - nu1_Hz) / ( log(vflux/bflux) - 3. * log(nu1_Hz/nu2_Hz) ); + if (verbose > 1) cout << "Tprime = " << tprime/1.38e-23*6.62e-34 << "\n"; + + + if( umag < -100. ){ // umag could not be read + if (verbose) cout << "Warning: star no. " << icatnum << + " has no Umag measurement. Calculating it from its Vmag = " << vmag << "\n"; + if (verbose) cout << " and Bmag = " << bmag << " assuming standard colour-colour-plot ... "; + + if((bmag-vmag) > 1.4){ + umag = bmag * 0.9; + } + else{ + if((bmag-vmag) > 0.5){ + umag = -0.5 + 1.37 * (bmag - vmag) + bmag; + } + else{ + if((bmag-vmag) <= 0.){ + umag = 4.07 * (bmag - vmag) + bmag; + } + else{ + umag = 0.175 * (bmag - vmag) + bmag; + } + } + } + + if (verbose) cout << " result Umag = " << umag << "\n"; + + if( umag < 5.0 ){ + cout << "Warning: star no. " << icatnum << " is bright (Vmag =" << vmag << ", Bmag = " + << bmag << ")\n and has no Umag measurement. Estimated Umag is "<< umag <<"\n"; + } + + } + else{ // umag available + if (verbose > 1) { + cout << "Test: star no. " << icatnum << + " has Umag = " << umag <<". Calculating it from its Vmag = " << vmag << "\n"; + cout << " and Bmag = " << bmag << " assuming standard colour-colour-plot ...\n "; + + if((bmag-vmag) > 1.4){ + xmag = bmag * 0.9; + } + else{ + if((bmag-vmag) > 0.5){ + xmag = -0.5 + 1.37 * (bmag - vmag) + bmag; + } + else{ + if((bmag-vmag) <= 0.){ + xmag = 4.07 * (bmag - vmag) + bmag; + } + else{ + xmag = 0.175 * (bmag - vmag) + bmag; + } + } + } + if (verbose > 2) + cout << "TEST " << umag <<" "<< xmag << " " << temp << " " << bmag << " " << vmag << "\n"; + + cout << " result Umag = " << xmag << "\n"; + } + + } + + if( rmag < -100. ){ // rmag not present (it's not part of the catalog) + + rflux = vflux * pow((VLMIN_nm + VLMAX_nm)/(RLMIN_nm + RLMAX_nm),3.) * + (exp(nu1_Hz/tprime) - 1.) / + (exp(LIGHTSPEED_mps/((RLMIN_nm+RLMAX_nm)/2.*1e-9)/tprime) - 1.); + + rmag = (log10(rflux) + 22.42)/(-0.4); + + } + + return(tprime); +} + + +//---------------------------------------------------------------------------- +// @name mag_nphot +// +// @desc translates magnitudes in number of photons, using log(flux)= -0.4*m+22.42 +// +//---------------------------------------------------------------------------- + +int star::mag_nphot(int np[4], float inttime_s, float radius_m, int verbose) { + + float bflux, vflux, uflux, rflux; + float bintensity, vintensity, uintensity, rintensity; + float unu1_Hz, unu2_Hz, bnu1_Hz, bnu2_Hz, vnu1_Hz, vnu2_Hz, rnu1_Hz, rnu2_Hz; + + // The flux is given in Watt/m2*Hz. + + uflux = pow(10.,-0.4*umag-22.42); + bflux = pow(10.,-0.4*bmag-22.42); + vflux = pow(10.,-0.4*vmag-22.42); + rflux = pow(10.,-0.4*rmag-22.42); + + if (verbose) cout<<"MAGS "< +#include +#include +#include +#include +#include "jcmacros.h" + +// Speed of light in m/s and Planck's constant in joules/s definition + +#ifndef LIGHTSPEED_mps +#define LIGHTSPEED_mps 3.000E8 +#endif + +#ifndef PLANCK_si +#define PLANCK_si 6.626E-34 +#endif + +#ifndef PI +#define PI 3.1415927 +#endif + + +// Wavelengths range for U, B and V magnitude + +#ifndef ULMIN_nm +#define ULMIN_nm 290.0 // the true limit is approx. 310 nm, but we extend this +#endif // in order to cover the complete PMT sensitivity + +#ifndef ULMAX_nm +#define ULMAX_nm 390.0 +#endif + +#ifndef BLMIN_nm +#define BLMIN_nm 390.0 +#endif + +#ifndef BLMAX_nm +#define BLMAX_nm 500.0 +#endif + +#ifndef VLMIN_nm +#define VLMIN_nm 500.0 +#endif + +#ifndef VLMAX_nm +#define VLMAX_nm 590.0 +#endif + +#ifndef RLMIN_nm +#define RLMIN_nm 590.0 +#endif + +#ifndef RLMAX_nm +#define RLMAX_nm 800.0 +#endif + +#ifndef SKY2000LINELENGTH +#define SKY2000LINELENGTH 520 +#endif + + + +class star{ + +public: + + int icatnum; + float ra_h, dec_deg, umag, bmag, vmag, rmag; + float u, v; + float ra_rad, dec_rad; + int numphot; + + star(); // constructor + int readstar(FILE *fp, int verbose); // read one star from the catalog, return TRUE + int printstar(); // write the parameters of one star to stdout + float calcmissingmags(int verbose);// calculate the magnitudes for those wavebands in which no data + // is available assuming a black body and using the V and B mags + int mag_nphot(int nph[4], float triggate_s, float radius_m, int verbose); + // translate the star magnitude into number of photons +}; + +#endif Index: trunk/MagicSoft/Simulation/Detector/Starfield/starfield.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/starfield.cxx (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/starfield.cxx (revision 341) @@ -0,0 +1,290 @@ +///////////////////////////////////////////////////////////////////////// +// Starfield Generator +// +// (c) 2000 D. Petry +// +///////////////////////////////////////////////////////////////////////// + + +#include "starfield.h" + +#define PROGRAMID "$Id: starfield.cxx,v 1.1 2000-01-21 13:36:36 petry Exp $" + +///////////////////////////////////////////////////////////////////////// + +main(int argc, char **argv) +{ + + parameters pars; + star stars[iMAXSTARS]; + photon photons[iMAXPHOT]; + + ifstream in; + FILE *catfile; + char parfilename[160]; + char catfilename[160]; + + + int i, j, k, numstars, subnumstars, starnumber, totalnumphot, photinside; + int totalphotinside, idum; + int istart_ra_h, iend_ra_h; + int nph[4]; // numbers of photons in the four wavebands + float lmin_nm[4] = {ULMIN_nm, BLMIN_nm, VLMIN_nm, RLMIN_nm}; // wave band definitions + float lmax_nm[4] = {ULMAX_nm, BLMAX_nm, VLMAX_nm, RLMAX_nm}; + float theta_rad, costheta, phi_rad, randtime, lambda_nm, xdum_m, ydum_m; + float angdist; + + //welcome + + cout << "This is STARFIELD. (c) 2000 D. Petry\n"; + cout << PROGRAMID << "\n"; + + // check command line arguments + + if(argc == 1){ + sprintf(parfilename, "starfield.par"); + } + else{ // a filename was given + sprintf(parfilename, "%s", argv[1]); + } + + // read parameter file + + in.open(parfilename, ios::in); // open the parameter file + + if(!in){ + cout << "Failed to open " << parfilename << "\n" + << "Value of stream \"in\" was " << in << "\n"; + cout << "\nThere shoud be a parameter file "<< parfilename + <<" in the working directory with the following format:\n-----\n"; + pars.usage(&cout); + cout << "-----\nExiting.\n"; + exit(1); + } + + cout << "Opened " << parfilename << " for reading ...\n"; + + if( !(pars.readparameters(&in)) ){ // read not OK? + if(!in.eof()){ // was the error not due to EOF? + cout << "Error: rdstate = " << in.rdstate() << "\n"; + } + else{ + cout << "Error: premature EOF in parameter file.\n"; + } + exit(1); + } + + in.close(); + + + // prepare loop over star catalog files + + cout << "SKY2000 - Master Star Catalog - Star Catalog Database, Version 2\n"; + cout << "Sande C.B., Warren W.H.Jr., Tracewell D.A., Home A.T., Miller A.C.\n"; + cout << "\n"; + + angdist = fmod( pars.catalog_fov_deg / cos( pars.ct_dec_rad ), 360.); + + if(angdist > 180.){ // too near to the pole, have to loop over all files + + istart_ra_h = 0; + iend_ra_h = 23; + + } + else{ // can loop over selected files only + + istart_ra_h = (int) (pars.ct_ra_h - angdist / 360. * 24.) - 1; + iend_ra_h = (int) (pars.ct_ra_h + angdist / 360. * 24. ) + 1; + + } + + //read catalog + + i = 0; + for (j = istart_ra_h; j <= iend_ra_h; j++){ + + subnumstars = 0; + + if ( j < 0){ + idum = j + 24; + } + else if ( j > 23 ){ + idum = j - 24; + } + else { + idum = j; + } + + sprintf(catfilename, "%s/sky%02d.dat", pars.datapath, idum); + + if((catfile = fopen(catfilename, "r")) == NULL){ // open the star catalog + cout << "Failed to open " << catfilename << "\n"; + exit(1); + } + cout << "Opened file " << catfilename << " for reading ...\n"; + + while( stars[i].readstar(catfile, pars.verbose) ){ // read next star OK + + angdist = acos( cos( pars.ct_dec_rad ) * cos( stars[i].dec_rad ) * + cos( stars[i].ra_rad - pars.ct_ra_rad ) + + sin( pars.ct_dec_rad ) * sin( stars[i].dec_rad ) ); + + if( angdist < pars.catalog_fov_deg / 180. * PI ){ // star in Field Of View? + + stars[i].calcmissingmags(pars.verbose); + if (pars.verbose) stars[i].printstar(); + + if( stars[i].umag > -100. ){ + i++; // accept star + subnumstars++; + } + else{ + cout << "Star rejected.\n"; + } + + if( i > iMAXSTARS ){ + i--; + cout << "Error: Star memory full. Accepted " << i << " stars.\n"; + break; + } + } + } + if( feof(catfile) ){ // was EOF reached? + cout << "EOF reached; accepted "<< subnumstars << " stars from this segment.\n"; + } + if( ferror(catfile) ){ // did an error occur? + cout << "Error while reading catalog file.\n"; + exit(1); + } + fclose(catfile); + if(i == iMAXSTARS){ + break; + } + } + + cout << "Accepted "<< i << " stars in total.\n"; + numstars = i; + + // loop over all photons from all stars, filling their fields + + totalnumphot=0; + + totalphotinside=0; + + for(i=0; i= 0. ){ + phi_rad = acos(stars[i].u); + } + else{ + phi_rad = 2.*PI - acos(stars[i].u); + } + + // calculate number of photons + + // mag_nphot() translates the star magnitude into number of photons, + // using the expression log(flux)=-0.4*m-22.42 for each waveband + // the resulting numbers ar stored in the array nph + + stars[i].mag_nphot(nph, pars.integtime_s, pars.mirr_radius_m, pars.verbose); + + // loop over all photons + + photinside=0; + + for(k=0; k < 4; k++){ // loop over wavebands + + for(j=0; j 2) + cout << "PH " << starnumber << " " << randtime << " " << xdum_m << " " + << ydum_m << " " << stars[i].u << " " << stars[i].v << " " << lambda_nm + << "\n"; + + photinside++; + + totalphotinside++; + + if(totalphotinside >= iMAXPHOT){ + cout << "Error: photon memory full. Can only store " << iMAXPHOT << " photons.\n"; + exit(1); + } + + totalnumphot++; + + } + else{ + + totalnumphot++; + continue; + + } // end if + + } // end photon loop + + } // end waveband loop + + stars[i].numphot = photinside; + + if (pars.verbose) cout<<"Star number= "< +#include +#include +#include +#include +#include +#include +#include + +// Include header files from the project. + +#include "convertcorsika.h" + +//#include "rand_time.h" + +//#include "rand_coord.h" + +#include "star.hxx" + +#include "photon.hxx" + +#include "jcmacros.h" + +#include "parameters.h" + +// Global constants. + + +// Maximum number of stars and maximum number of photons for +// memory allocation + + +#ifndef iMAXSTARS +#define iMAXSTARS 5000 +#endif + +#ifndef iMAXPHOT +#define iMAXPHOT 250000 +#endif + +#ifndef PI +#define PI 3.1415927 +#endif + + +//Seed for the random generators. + +long int idum=-1; + + +//Function definitions + + +//Random generator functions. + +float rand_un_gen(long *idum); + + + +inline float rand_coord(float radius){ +// generate a coordinate on the floor for every photon arriving from the star. + + float p_m; + float randc; + + randc=rand_un_gen(&idum); + + p_m= 2.0*radius * randc - radius; + + return (p_m); +} + +// generate for every photon a random arrival time within a time window +// defined by the intgreation time. + +inline float rand_time(float inttime_s){ + + float randt; + + randt= inttime_s * rand_un_gen(&idum); + + return(randt); +} + +// The rand_lambda function will generate for every photon a random wavelength +// in the given waveband + +inline float rand_lambda(float lmin, float lmax){ + + float randt; + + randt=(lmax - lmin) * rand_un_gen(&idum) + lmin ; + + return(randt); +} + + +#endif + Index: trunk/MagicSoft/Simulation/Detector/Starfield/starfield.par =================================================================== --- trunk/MagicSoft/Simulation/Detector/Starfield/starfield.par (revision 341) +++ trunk/MagicSoft/Simulation/Detector/Starfield/starfield.par (revision 341) @@ -0,0 +1,13 @@ +Starfield Generator Parameters, Date: 21-1-2000, Comment: Example (Crab Nebula) +Center of FOV: ira_hours ira_min ira_sec idec_degrees idec_arcmin dec_arcsec: +05 34 32 +22 00 52.1 +Radius of the FOV for the catalog readout (degrees): +2.0 +Integration time for the calculation of the number of photons (seconds): +50e-9 +Mirror radius for the generation of random impact points (meters): +10.0 +Path inside which the star catalog data can be found: +/usr/users/petry/temp/starfield/sky2000 +Verbosity level (0 = not verbose, 1 = verbose, 2 = very verbose, 3 = very very ...): +0