Index: trunk/MagicSoft/Simulation/Detector/Reflector/Makefile =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/Makefile (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/Makefile (revision 292) @@ -0,0 +1,251 @@ +################################################################## +# +# makefile +# +# @file makefile +# @title Simulation of the camera and trigger logic +# @author J C Gonz\'alez +# @email gonzalez@mppmu.mpg.de +# @date Fri Mar 12 11:51:11 MET 1999 +# +#_______________________________________________________________ +# +# Created: Fri Mar 12 11:51:11 MET 1999 +# Author: Jose Carlos Gonzalez +# Purpose: Makefile for the compilation of the camera program +# Notes: +# +#--------------------------------------------------------------- +# +# $RCSfile: Makefile,v $ +# $Revision: 1.1.1.1 $ +# $Author: harald $ +# $Date: 1999-10-29 07:00:32 $ +# +################################################################## +# @maintitle + +# @code + +INCLUDEMK = config.mk +include ${INCLUDEMK} + +# @endcode + +# @code + +# common flags +INCLUDES = -I${INCLUDE} \ + -I${INCLUDE_COR} \ + -I${INCLUDE_MC} \ + -I${INCLUDE_REFL} \ + -I/usr/include + +#CERNLIBDIR = ${CERNDIR}/pro/lib/ +#CERNLIB = -L${CERNLIBDIR} -lgraflib -lgrafX11 -lpacklib -lkernlib -lpawlib + +RANLIB = -L${RANLIBDIR} -lranlib + +# special flags + +osf_FORLIBS = -lUfor -lfor -lutil -lots -lm +linux_FORLIBS = -lm +generic_FORLIBS = -lm + +FORLIBS = ${${SYSTEM}_FORLIBS} + +# compilation and linking flags + +CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} ${DEBUG} +CFLAGS = ${CXXFLAGS} +FFLAGS = ${CXXFLAGS} +LIBS = ${RANLIB} ${FORLIBS} +#LIBS = ${CERNLIB} ${RANLIB} ${FORLIBS} + +#------------------------------------------------------------------------------ + +#.SILENT: + +.SUFFIXES: .c .cxx .C .c++ .h .hxx .H .h++ .o .so .f + +SRCS = \ + ${INCLUDE_COR}/COREventHeader.cxx \ + ${INCLUDE_COR}/CORParticle.cxx \ + ${INCLUDE_COR}/CORStatfile.cxx \ + ${INCLUDE_MC}/MCEventHeader.cxx \ + ${INCLUDE_MC}/MCCphoton.cxx \ + atm.cxx \ + attenu.f \ + readparam.cxx \ + reflector.cxx + +HEADERS = \ + COREventHeader.hxx \ + CORParticle.hxx \ + CORStatfile.hxx \ + MCEventHeader.hxx \ + MCCphoton.hxx \ + lagrange.h \ + atm.h \ + readparam.h \ + reflector.h + +OBJS = \ + ${INCLUDE_COR}/COREventHeader.o \ + ${INCLUDE_COR}/CORParticle.o \ + ${INCLUDE_COR}/CORStatfile.o \ + ${INCLUDE_MC}/MCEventHeader.o \ + ${INCLUDE_MC}/MCCphoton.o \ + atm.o \ + attenu.o \ + readparam.o \ + reflector.o + +PROGRAM=reflector + +############################################################ + +all: ${PROGRAM} + +depend: + @makedepend $(SRCS) -fMakefile 2> /dev/null + +doc: reflector-doc + +reflector-doc: + @echo "Generating documentation for camera . . . " + $(DOCUM) -latex -o reflector.tex \ + reflector.cxx reflector.h \ + readparam.cxx readparam.h \ + atm.cxx atm.h + latex "\nonstopmode\input{reflector.tex}" && \ + makeindex reflector && \ + latex "\nonstopmode\input{reflector.tex}" && \ + latex "\nonstopmode\input{reflector.tex}" + @echo "Files reflector.tex and reflector.dvi generated." + +rate: + @echo "Rating documentation inside code . . . " + $(RATE) \ + reflector.cxx reflector.h \ + readparam.cxx readparam.h \ + atm.cxx atm.h + +${PROGRAM}: $(OBJS) + @echo "Linking..." $@ + $(CXX) $(CXXFLAGS) $(OBJS) $(LIBS) -o $@ + @echo "done." + +.cxx.o: + @echo "Compiling " $< + $(CXX) $(CXXFLAGS) -c $< -o $@ + +.c.o: + @echo "Compiling " $< + $(CC) $(CFLAGS) -c $< -o $@ + +.f.o: + @echo "Compiling " $< + $(F77) $(FFLAGS) -c $< -o $@ + +lclean: + @echo "Cleanning..." + @rm -f *.o core + +clean: + @echo "Cleanning..." + @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 . . ." + @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. + +../include-CORSIKA/COREventHeader.o: ../include-CORSIKA/COREventHeader.hxx +../include-CORSIKA/COREventHeader.o: /usr/include/stdlib.h +../include-CORSIKA/COREventHeader.o: /usr/include/standards.h +../include-CORSIKA/COREventHeader.o: /usr/include/getopt.h +../include-CORSIKA/COREventHeader.o: /usr/include/sys/types.h +../include-CORSIKA/COREventHeader.o: /usr/include/mach/machine/vm_types.h +../include-CORSIKA/COREventHeader.o: /usr/include/sys/select.h +../include-CORSIKA/COREventHeader.o: /usr/include/math.h +../include-CORSIKA/CORParticle.o: ../include-CORSIKA/CORParticle.hxx +../include-CORSIKA/CORParticle.o: /usr/include/stdlib.h +../include-CORSIKA/CORParticle.o: /usr/include/standards.h +../include-CORSIKA/CORParticle.o: /usr/include/getopt.h +../include-CORSIKA/CORParticle.o: /usr/include/sys/types.h +../include-CORSIKA/CORParticle.o: /usr/include/mach/machine/vm_types.h +../include-CORSIKA/CORParticle.o: /usr/include/sys/select.h +../include-CORSIKA/CORParticle.o: /usr/include/math.h +../include-CORSIKA/CORStatfile.o: ../include-CORSIKA/CORStatfile.hxx +../include-CORSIKA/CORStatfile.o: /usr/include/stdlib.h +../include-CORSIKA/CORStatfile.o: /usr/include/standards.h +../include-CORSIKA/CORStatfile.o: /usr/include/getopt.h +../include-CORSIKA/CORStatfile.o: /usr/include/sys/types.h +../include-CORSIKA/CORStatfile.o: /usr/include/mach/machine/vm_types.h +../include-CORSIKA/CORStatfile.o: /usr/include/sys/select.h +../include-CORSIKA/CORStatfile.o: /usr/include/math.h +../include-MC/MCEventHeader.o: ../include-MC/MCEventHeader.hxx +../include-MC/MCEventHeader.o: /usr/include/stdlib.h /usr/include/standards.h +../include-MC/MCEventHeader.o: /usr/include/getopt.h /usr/include/sys/types.h +../include-MC/MCEventHeader.o: /usr/include/mach/machine/vm_types.h +../include-MC/MCEventHeader.o: /usr/include/sys/select.h /usr/include/math.h +../include-MC/MCEventHeader.o: ../include-CORSIKA/COREventHeader.hxx +../include-MC/MCCphoton.o: ../include-MC/MCCphoton.hxx /usr/include/stdlib.h +../include-MC/MCCphoton.o: /usr/include/standards.h /usr/include/getopt.h +../include-MC/MCCphoton.o: /usr/include/sys/types.h +../include-MC/MCCphoton.o: /usr/include/mach/machine/vm_types.h +../include-MC/MCCphoton.o: /usr/include/sys/select.h /usr/include/string.h +../include-MC/MCCphoton.o: /usr/include/strings.h /usr/include/math.h +atm.o: atm.h /usr/include/stdlib.h /usr/include/standards.h +atm.o: /usr/include/getopt.h /usr/include/sys/types.h +atm.o: /usr/include/mach/machine/vm_types.h /usr/include/sys/select.h +atm.o: /usr/include/stdio.h /usr/include/sys/seek.h /usr/include/va_list.h +atm.o: /usr/include/sys/limits.h /usr/include/sys/machine/machlimits.h +atm.o: /usr/include/sys/syslimits.h /usr/include/sys/machine/machtime.h +atm.o: /usr/include/sys/rt_limits.h /usr/include/string.h +atm.o: /usr/include/strings.h /usr/include/math.h reflector-v.h +readparam.o: readparam.h /usr/include/stdlib.h /usr/include/standards.h +readparam.o: /usr/include/getopt.h /usr/include/sys/types.h +readparam.o: /usr/include/mach/machine/vm_types.h /usr/include/sys/select.h +readparam.o: /usr/include/stdio.h /usr/include/sys/seek.h +readparam.o: /usr/include/va_list.h /usr/include/sys/limits.h +readparam.o: /usr/include/sys/machine/machlimits.h +readparam.o: /usr/include/sys/syslimits.h /usr/include/sys/machine/machtime.h +readparam.o: /usr/include/sys/rt_limits.h /usr/include/string.h +readparam.o: /usr/include/strings.h /usr/include/math.h /usr/include/float.h +readparam.o: /usr/include/fp_class.h /usr/include/unistd.h +readparam.o: /usr/include/sys/access.h reflector-v.h atm.h +reflector.o: reflector.h /usr/include/stdlib.h /usr/include/standards.h +reflector.o: /usr/include/getopt.h /usr/include/sys/types.h +reflector.o: /usr/include/mach/machine/vm_types.h /usr/include/sys/select.h +reflector.o: /usr/include/stdio.h /usr/include/sys/seek.h +reflector.o: /usr/include/va_list.h /usr/include/sys/limits.h +reflector.o: /usr/include/sys/machine/machlimits.h +reflector.o: /usr/include/sys/syslimits.h /usr/include/sys/machine/machtime.h +reflector.o: /usr/include/sys/rt_limits.h /usr/include/string.h +reflector.o: /usr/include/strings.h /usr/include/stdarg.h /usr/include/math.h +reflector.o: /usr/include/dirent.h /usr/include/unistd.h +reflector.o: /usr/include/sys/access.h reflector-v.h readparam.h +reflector.o: /usr/include/float.h /usr/include/fp_class.h atm.h +reflector.o: ../include-CORSIKA/COREventHeader.hxx +reflector.o: ../include-CORSIKA/CORParticle.hxx +reflector.o: ../include-CORSIKA/CORStatfile.hxx +reflector.o: ../include-MC/MCEventHeader.hxx ../include-MC/MCCphoton.hxx +reflector.o: /usr/include/ranlib.h /usr/include/ar.h Index: trunk/MagicSoft/Simulation/Detector/Reflector/atm.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/atm.cxx (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/atm.cxx (revision 292) @@ -0,0 +1,225 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= atm +//= +//= @file atm.cxx +//= @desc Header file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Simulation of atmospheric absorption +//= Notes: +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: atm.cxx,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-10-29 07:00:32 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |atm.cxx|. + +/*!@" + + This section describes briefly the source code for the file + |atm.cxx|. + + @"*/ + +//!@subsection Include files. + +/*!@" + + All the defines are located in the file {\tt atm.h}. + + @"*/ + +//!@{ + +#include "atm.h" + +//!@} + +//!@subsection Definition of global variables. + +//!@{ + +static ATM_MODELS_TYPE Atm_mod_N = ATM_NOATMOSPHERE; +static char Atm_mod[ATM_NAME_MAX_LENGTH]; + +//!@} + +//!@subsection The function |atm()|. + +//!--------------------------------------------------------------------- +// @name atm +// +// @desc simulation of atmospheric absorption +// +// @var wavelength Wavelength of the photon (nm) +// @var height Height of generation (cm) +// @var theta Zenith angle of trajectory(rad) +// +// @return Transmission coefficient +// +// @date Mon Sep 14 13:27:56 MET DST 1997 +//---------------------------------------------------------------------- +// @function + +//!@{ +float +atm(float wavelength, float height, float theta) +{ + + float transmittance; //@< final transmittance of the atm. (returned) + + transmittance = 1.0; + + switch ( Atm_mod_N ) { + + case ATM_NOATMOSPHERE: //@< no atmosphere at all: transmittance = 100% + + break; + + case ATM_90PERCENT: //@< atmosphere with transmittance = 90% + + transmittance = 0.9; + break; + + case ATM_ISOTHERMAL: //@< isothermal approximation + + // here the transmittance function + + break; + + case ATM_CORSIKA: //@< atmosphere as defined in CORSIKA + + attenu( &wavelength, &height, &theta, &transmittance); + + break; + + } + + return ( transmittance ); + +} +//!@} + +//!--------------------------------------------------------------------- +// @name set_atm_mod +// +// @desc set model of atmospheric absorption +// +// @var i Model of atmospheric absorption to be used +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +void +set_atm_mod(ATM_MODELS_TYPE i) +{ + Atm_mod_N = i; + strcpy(Atm_mod, ATM_MODELS_NAMES[ Atm_mod_N ]); +} +//!@} + +//!--------------------------------------------------------------------- +// @name conv_atm_mod +// +// @desc convert string with name of atm.model to number +// +// @var *s Name of model of atmospheric absorption +// +// @return Number for that model +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +ATM_MODELS_TYPE +conv_atm_mod(char *s) +{ + int i; + char s2[40]; + + // look at each model + for (i=0; i<=ATM_CORSIKA; i++) { + strcpy(s2, ATM_MODELS_NAMES[i]); + if (strcmp(s, s2) == NULL) + break; + } + + // return the name of the model + return ( (ATM_MODELS_TYPE) ((i<=ATM_CORSIKA) ? i : -1) ); + + // return ( (ATM_MODELS_TYPE) 0 ); +} +//!@} + +//!--------------------------------------------------------------------- +// @name get_atm_mod +// +// @desc get name of the atmospheric model to be used +// +// @return Name of the current atmospheric model to be used +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +char * +get_atm_mod(void) +{ + return (Atm_mod); +} +//!@} + +//!--------------------------------------------------------------------- +// @name get_atm_mod_n +// +// @desc get name of the atmospheric model number i +// +// @var i Number of atmospheric model +// +// @return Name of corresponding atmospheric model +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +char * +get_atm_mod_n(ATM_MODELS_TYPE i) +{ + return( ATM_MODELS_NAMES[ Atm_mod_N ] ); +} +//!@} + + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ +// +// $Log: not supported by cvs2svn $ +// Revision 1.3 1999/03/24 16:32:59 gonzalez +// REFLECTOR 1.1: Release +// +// +// +//!@} + +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Reflector/atm.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/atm.h (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/atm.h (revision 292) @@ -0,0 +1,145 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= atm +//= +//= @file atm.cxx +//= @desc Header file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Simulation of atmospheric absorption +//= Notes: +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: atm.h,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-10-29 07:00:32 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |atm.h|. + +/*!@" + + In this section you can find the source code for the file |atm.h|. + + @"*/ + +//!@{ +#ifndef _atm_ +#define _atm_ + +#ifdef _this_ +#define _define_extern_ +#endif + +#ifndef _this_ +#define _this_ atm +#endif +//!@} + +//!@subsection Include files. + +//!@{ + +#include +#include +#include +#include +#include +#include + +#include "jcmacros.h" +#include "jcdebug.h" + +#include "reflector-v.h" + +//!@} + +//!@subsection Macro-definitions, and constants. + +//!@{ + +// now we define the list ITEM_LIST of possible items in +// the parameters file. note that they are functions of +// T(x). we will change T(x) to be the real item or the +// string of this item when needed + +#define ATM_MODELS_LIST /* LIST OF POSIBLE MODELS OF ATMOSPHERE */ \ +T(ATM_NOATMOSPHERE), /* no atmosphere at all: transmittance = 100% */ \ +T(ATM_90PERCENT), /* atmosphere with transmittance = 90% */ \ +T(ATM_ISOTHERMAL), /* isothermal approximation */ \ +T(ATM_CORSIKA) /* atmosphere as defined in CORSIKA */ + +#define T(x) x // define T() as the name as it is + +enum ATM_MODELS_TYPE { + ATM_MODELS_LIST +}; + +#undef T + +#define T(x) #x // define T() as the string of x + +char *const ATM_MODELS_NAMES[] = { + ATM_MODELS_LIST +}; + +#undef T + +#define ATM_NAME_MAX_LENGTH 40 + +float atm(float wavelength, float height, float theta); +#ifdef _define_extern_ +extern void set_atm_mod(ATM_MODELS_TYPE i); +extern ATM_MODELS_TYPE conv_atm_mod(char *atmod); +#else +void set_atm_mod(ATM_MODELS_TYPE i); +ATM_MODELS_TYPE conv_atm_mod(char *atmod); +#endif +char *get_atm_mod(); +char *get_atm_mod_n(ATM_MODELS_TYPE i); + + +// routines written in Fortran + +#define attenu attenu_ +extern "C" { + void attenu(float *, float *, float *, float *); +} + +//!@} + +//!@{ + +#endif // ! _atm_ + +//!@} + + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ + +/* + * $Log: not supported by cvs2svn $ + * Revision 1.3 1999/03/24 16:33:02 gonzalez + * REFLECTOR 1.1: Release + * + * Revision 1.2 1999/01/21 16:03:32 gonzalez + * Only small modifications + * + */ + +//!@} +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Reflector/attenu.f =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/attenu.f (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/attenu.f (revision 292) @@ -0,0 +1,422 @@ +********************************************************************* +* * +* Atmospheric Absorption for Rayleigh, Mie and Ozone. * +* * +* Created: May-98 * +* Author: Aitor Ibarra Ibaibarriaga * +* Jose Carlos Gonzalez * +* * +********************************************************************* + +c @T \newpage + +c @section Source code of {\tt attenu.f} + +* @text +* Copyright $\copyright$ 1998, Aitor Ibarra Ibaibarriaga +* @endtext +c @code + + SUBROUTINE attenu(wavelength, height, theta, tr_atmos) + +C----------------------------------------------------------------------C +C RHO (DENSITY) F(UNCTION) C +C C +C CALCULATES DENSITY (G/CM**3) OF ATMOSPHERE DEPENDING ON HEIGHT (CM) C +C (US STANDARD ATMOSPHERE) C +C THIS FUNCTION IS CALLED FROM ININKG, UPDATE, CERENE, CERENH C +C ARGUMENT: C +C ARG = HEIGHT IN CM C +C----------------------------------------------------------------------C + + COMMON /ATMOS/ AATM,BATM,CATM,LAHG,RHOF,LONG,OZ_ABI,AE_ABI + DOUBLE PRECISION AATM(5),BATM(5),CATM(5),LAHG(5),RHOF(5),Lm(12) + DOUBLE PRECISION H,RT, m,OZ_ABI(48,12),AE_ABI(28,12) + DOUBLE PRECISION XR, TOT_OZ, TOT_AE, T_Ray,T_Mie, T_Oz, + + RHO_TOT, RHO_FI, RHOFP + + REAL wavelength, height, theta, tr_atmos + real trr, trm, tro + REAL LONG(12) + INTEGER I,CON_OZ,CON_MI J, ROW + + DATA OZ_ABI / + + 0.2880000D+00,0.5405000D+00,0.7775000D+00,0.1009000D+01, + + 0.1241500D+01,0.1480500D+01,0.1750500D+01,0.2085000D+01, + + 0.2514500D+01,0.3087500D+01,0.3864500D+01,0.4817500D+01, + + 0.5847500D+01,0.6917500D+01,0.8052500D+01,0.9287499D+01, + + 0.1068750D+02,0.1231250D+02,0.1415750D+02,0.1617750D+02, + + 0.1827250D+02,0.2034750D+02,0.2232750D+02,0.2414750D+02, + + 0.2575750D+02,0.2715250D+02,0.2836750D+02,0.2941100D+02, + + 0.3031000D+02,0.3109250D+02,0.3176300D+02,0.3232750D+02, + + 0.3281200D+02,0.3323200D+02,0.3358350D+02,0.3387750D+02, + + 0.3412650D+02,0.3434000D+02,0.3451900D+02,0.3466250D+02, + + 0.3477480D+02,0.3486355D+02,0.3493355D+02,0.3498775D+02, + + 0.3503010D+02,0.3506360D+02,0.3509020D+02,0.3511185D+02, + + 0.2740000D-01,0.5140000D-01,0.7395000D-01,0.9600000D-01, + + 0.1181500D+00,0.1409000D+00,0.1666000D+00,0.1984500D+00, + + 0.2393500D+00,0.2939500D+00,0.3679500D+00,0.4589500D+00, + + 0.5573000D+00,0.6593000D+00,0.7673000D+00,0.8848000D+00, + + 0.1017800D+01,0.1172300D+01,0.1348300D+01,0.1540800D+01, + + 0.1740300D+01,0.1937800D+01,0.2126300D+01,0.2299800D+01, + + 0.2453300D+01,0.2586300D+01,0.2702300D+01,0.2801900D+01, + + 0.2887550D+01,0.2962050D+01,0.3025900D+01,0.3079700D+01, + + 0.3125850D+01,0.3165850D+01,0.3199350D+01,0.3227400D+01, + + 0.3251150D+01,0.3271500D+01,0.3288600D+01,0.3302300D+01, + + 0.3312995D+01,0.3321445D+01,0.3328110D+01,0.3333270D+01, + + 0.3337305D+01,0.3340500D+01,0.3343035D+01,0.1334510D+01, + + 0.2435000D-02,0.4570000D-02,0.6575000D-02,0.8535000D-02, + + 0.1050500D-01,0.1253000D-01,0.1481500D-01,0.1764500D-01, + + 0.2128000D-01,0.2613500D-01,0.3272000D-01,0.4081000D-01, + + 0.4957000D-01,0.5866000D-01,0.6827000D-01,0.7875500D-01, + + 0.9065500D-01,0.1044050D+00,0.1200050D+00,0.1371050D+00, + + 0.1548550D+00,0.1724050D+00,0.1891550D+00,0.2045550D+00, + + 0.2182050D+00,0.2300550D+00,0.2403600D+00,0.2492200D+00, + + 0.2568350D+00,0.2634550D+00,0.2691300D+00,0.2739150D+00, + + 0.2780200D+00,0.2815750D+00,0.2845500D+00,0.2870400D+00, + + 0.2891500D+00,0.2909600D+00,0.2924750D+00,0.2936900D+00, + + 0.2946425D+00,0.2953940D+00,0.2959865D+00,0.2964455D+00, + + 0.2968045D+00,0.2970885D+00,0.2973140D+00,0.2974975D+00, + + 0.1740000D-03,0.3265000D-03,0.4695000D-03,0.6090000D-03, + + 0.7495000D-03,0.8940000D-03,0.1057000D-02,0.1259000D-02, + + 0.1518000D-02,0.1863500D-02,0.2332500D-02,0.2909000D-02, + + 0.3533000D-02,0.4180500D-02,0.4865000D-02,0.5610500D-02, + + 0.6455500D-02,0.7435000D-02,0.8550000D-02,0.9770000D-02, + + 0.1103500D-01,0.1229000D-01,0.1348500D-01,0.1458000D-01, + + 0.1555100D-01,0.1639550D-01,0.1713150D-01,0.1776300D-01, + + 0.1830600D-01,0.1877800D-01,0.1918200D-01,0.1952250D-01, + + 0.1981500D-01,0.2006850D-01,0.2028050D-01,0.2045801D-01, + + 0.2060851D-01,0.2073750D-01,0.2084566D-01,0.2093241D-01, + + 0.2100026D-01,0.2105381D-01,0.2109606D-01,0.2112876D-01, + + 0.2115431D-01,0.2117456D-01,0.2119066D-01,0.2120376D-01, + + 0.4885000D-05,0.9170000D-05,0.1319500D-04,0.1713000D-04, + + 0.2108000D-04,0.2513500D-04,0.2971500D-04,0.3539500D-04, + + 0.4268500D-04,0.5242500D-04,0.6562500D-04,0.8182500D-04, + + 0.9937500D-04,0.1175750D-03,0.1368250D-03,0.1578250D-03, + + 0.1816250D-03,0.2091750D-03,0.2404750D-03,0.2747750D-03, + + 0.3103250D-03,0.3454750D-03,0.3790250D-03,0.4098750D-03, + + 0.4372250D-03,0.4609750D-03,0.4816750D-03,0.4994750D-03, + + 0.5147750D-03,0.5280750D-03,0.5394750D-03,0.5490700D-03, + + 0.5572950D-03,0.5644250D-03,0.5703950D-03,0.5753899D-03, + + 0.5796200D-03,0.5832500D-03,0.5862950D-03,0.5887350D-03, + + 0.5906400D-03,0.5921450D-03,0.5933350D-03,0.5942565D-03, + + 0.5949755D-03,0.5955440D-03,0.5959955D-03,0.5963635D-03, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00, + + 0.9525000D-05,0.1785500D-04,0.2567500D-04,0.3332000D-04, + + 0.4100000D-04,0.4889000D-04,0.5779500D-04,0.6881000D-04, + + 0.8296000D-04,0.1018600D-03,0.1275100D-03,0.1590600D-03, + + 0.1932100D-03,0.2286100D-03,0.2660100D-03,0.3067600D-03, + + 0.3529600D-03,0.4065100D-03,0.4674100D-03,0.5341100D-03, + + 0.6032600D-03,0.6717100D-03,0.7370100D-03,0.7970100D-03, + + 0.8501600D-03,0.8963600D-03,0.9366100D-03,0.9711101D-03, + + 0.1000810D-02,0.1026660D-02,0.1048810D-02,0.1067460D-02, + + 0.1083460D-02,0.1097310D-02,0.1108910D-02,0.1118640D-02, + + 0.1126865D-02,0.1133915D-02,0.1139830D-02,0.1144570D-02, + + 0.1148275D-02,0.1151200D-02,0.1153510D-02,0.1155300D-02, + + 0.1156700D-02,0.1162205D-02,0.1170990D-02,0.1178145D-02, + + 0.9360000D-04,0.1757000D-03,0.2528000D-03,0.3281500D-03, + + 0.4038500D-03,0.4816500D-03,0.5694500D-03,0.6784000D-03, + + 0.8184000D-03,0.1004900D-02,0.1257900D-02,0.1568900D-02, + + 0.1905400D-02,0.2254400D-02,0.2623400D-02,0.3025400D-02, + + 0.3480900D-02,0.4008900D-02,0.4609400D-02,0.5266900D-02, + + 0.5948400D-02,0.6622900D-02,0.7266400D-02,0.7857900D-02, + + 0.8381900D-02,0.8836901D-02,0.9233401D-02,0.9573901D-02, + + 0.9866901D-02,0.1012140D-01,0.1033940D-01,0.1052340D-01, + + 0.1068140D-01,0.1081840D-01,0.1093290D-01,0.1102855D-01, + + 0.1110965D-01,0.1117920D-01,0.1123750D-01,0.1128425D-01, + + 0.1132085D-01,0.1134975D-01,0.1137255D-01,0.1139020D-01, + + 0.1140400D-01,0.1141492D-01,0.1142358D-01,0.1143063D-01, + + 0.2500000D-03,0.4690000D-03,0.6745000D-03,0.8755000D-03, + + 0.1077500D-02,0.1285000D-02,0.1519500D-02,0.1810000D-02, + + 0.2182500D-02,0.2679500D-02,0.3353500D-02,0.4182000D-02, + + 0.5079000D-02,0.6010000D-02,0.6994000D-02,0.8064000D-02, + + 0.9279000D-02,0.1068900D-01,0.1228900D-01,0.1403900D-01, + + 0.1585400D-01,0.1765400D-01,0.1937400D-01,0.2095400D-01, + + 0.2235400D-01,0.2356900D-01,0.2462600D-01,0.2553350D-01, + + 0.2631400D-01,0.2699250D-01,0.2757350D-01,0.2806300D-01, + + 0.2848350D-01,0.2884800D-01,0.2915300D-01,0.2940850D-01, + + 0.2962500D-01,0.2981050D-01,0.2996600D-01,0.3009050D-01, + + 0.3018780D-01,0.3026480D-01,0.3032550D-01,0.3037250D-01, + + 0.3040925D-01,0.3043835D-01,0.3046145D-01,0.3048025D-01, + + 0.3585000D-03,0.6725000D-03,0.9675000D-03,0.1256000D-02, + + 0.1545500D-02,0.1843000D-02,0.2179000D-02,0.2595500D-02, + + 0.3130000D-02,0.3843500D-02,0.4813500D-02,0.6003500D-02, + + 0.7288500D-02,0.8623499D-02,0.1003850D-01,0.1157850D-01, + + 0.1331850D-01,0.1533350D-01,0.1762850D-01,0.2014350D-01, + + 0.2274850D-01,0.2532850D-01,0.2779350D-01,0.3005850D-01, + + 0.3206350D-01,0.3380350D-01,0.3531851D-01,0.3661850D-01, + + 0.3773851D-01,0.3871351D-01,0.3954751D-01,0.4025051D-01, + + 0.4085401D-01,0.4137701D-01,0.4181501D-01,0.4218151D-01, + + 0.4249151D-01,0.4275751D-01,0.4298101D-01,0.4316001D-01, + + 0.4330001D-01,0.4341061D-01,0.4349771D-01,0.4356516D-01, + + 0.4361791D-01,0.4365961D-01,0.4369271D-01,0.4371971D-01, + + 0.1685000D-03,0.3160000D-03,0.4545000D-03,0.5900000D-03, + + 0.7260000D-03,0.8655000D-03,0.1023000D-02,0.1218500D-02, + + 0.1469500D-02,0.1804500D-02,0.2259000D-02,0.2817500D-02, + + 0.3422000D-02,0.4049500D-02,0.4712999D-02,0.5434999D-02, + + 0.6252999D-02,0.7203000D-02,0.8283000D-02,0.9462999D-02, + + 0.1068800D-01,0.1190300D-01,0.1306300D-01,0.1412800D-01, + + 0.1507000D-01,0.1588850D-01,0.1660150D-01,0.1721300D-01, + + 0.1773900D-01,0.1819650D-01,0.1858850D-01,0.1891850D-01, + + 0.1920150D-01,0.1944700D-01,0.1965250D-01,0.1982450D-01, + + 0.1997050D-01,0.2009550D-01,0.2020010D-01,0.2028410D-01, + + 0.2034985D-01,0.2040175D-01,0.2044265D-01,0.2047435D-01, + + 0.2049915D-01,0.2051875D-01,0.2053430D-01,0.2054695D-01 / + + DATA AE_ABI / + + 0.3645000E-01,0.5211000E-01,0.5913000E-01,0.6203000E-01, + + 0.6304000E-01,0.6340450E-01,0.6353275E-01,0.6358405E-01, + + 0.6361780E-01,0.6364885E-01,0.6367990E-01,0.6371365E-01, + + 0.6375955E-01,0.6383380E-01,0.6392965E-01,0.6403360E-01, + + 0.6414810E-01,0.6426660E-01,0.6438010E-01,0.6448960E-01, + + 0.6459510E-01,0.6468170E-01,0.6474110E-01,0.6478295E-01, + + 0.6481670E-01,0.6484775E-01,0.6487610E-01,0.6490240E-01, + + 0.3510000E-01,0.5018000E-01,0.5694000E-01,0.5973500E-01, + + 0.6071000E-01,0.6106100E-01,0.6118450E-01,0.6123390E-01, + + 0.6126640E-01,0.6129630E-01,0.6132620E-01,0.6135870E-01, + + 0.6140290E-01,0.6147440E-01,0.6156670E-01,0.6166680E-01, + + 0.6177730E-01,0.6189180E-01,0.6200130E-01,0.6210680E-01, + + 0.6220820E-01,0.6229140E-01,0.6234860E-01,0.6238890E-01, + + 0.6242140E-01,0.6245130E-01,0.6247860E-01,0.6250395E-01, + + 0.3375000E-01,0.4825000E-01,0.5475000E-01,0.5744000E-01, + + 0.5838000E-01,0.5871750E-01,0.5883625E-01,0.5888375E-01, + + 0.5891500E-01,0.5894375E-01,0.5897250E-01,0.5900375E-01, + + 0.5904625E-01,0.5911500E-01,0.5920375E-01,0.5930000E-01, + + 0.5940650E-01,0.5951700E-01,0.5962200E-01,0.5972300E-01, + + 0.5982050E-01,0.5990050E-01,0.5995550E-01,0.5999425E-01, + + 0.6002550E-01,0.6005425E-01,0.6008050E-01,0.6010490E-01, + + 0.3240000E-01,0.4632000E-01,0.5256000E-01,0.5514000E-01, + + 0.5604000E-01,0.5636400E-01,0.5647800E-01,0.5652360E-01, + + 0.5655360E-01,0.5658120E-01,0.5660880E-01,0.5663880E-01, + + 0.5667960E-01,0.5674560E-01,0.5683080E-01,0.5692320E-01, + + 0.5702520E-01,0.5713070E-01,0.5723140E-01,0.5732860E-01, + + 0.5742220E-01,0.5749900E-01,0.5755180E-01,0.5758900E-01, + + 0.5761900E-01,0.5764660E-01,0.5767180E-01,0.5769520E-01, + + 0.3240000E-01,0.4632000E-01,0.5256000E-01,0.5514000E-01, + + 0.5604000E-01,0.5636400E-01,0.5647800E-01,0.5652360E-01, + + 0.5655360E-01,0.5658120E-01,0.5660880E-01,0.5663880E-01, + + 0.5667960E-01,0.5674560E-01,0.5683080E-01,0.5692320E-01, + + 0.5702520E-01,0.5713070E-01,0.5723140E-01,0.5732860E-01, + + 0.5742220E-01,0.5749900E-01,0.5755180E-01,0.5758900E-01, + + 0.5761900E-01,0.5764660E-01,0.5767180E-01,0.5769520E-01, + + 0.3105000E-01,0.4439000E-01,0.5037000E-01,0.5284500E-01, + + 0.5371000E-01,0.5402050E-01,0.5412975E-01,0.5417345E-01, + + 0.5420220E-01,0.5422865E-01,0.5425510E-01,0.5428385E-01, + + 0.5432295E-01,0.5438620E-01,0.5446785E-01,0.5455640E-01, + + 0.5465390E-01,0.5475485E-01,0.5485145E-01,0.5494460E-01, + + 0.5503430E-01,0.5510790E-01,0.5515850E-01,0.5519415E-01, + + 0.5522290E-01,0.5524935E-01,0.5527350E-01,0.5529590E-01, + + 0.2700000E-01,0.3860000E-01,0.4380000E-01,0.4595000E-01, + + 0.4670000E-01,0.4697000E-01,0.4706500E-01,0.4710300E-01, + + 0.4712800E-01,0.4715100E-01,0.4717400E-01,0.4719900E-01, + + 0.4723300E-01,0.4728800E-01,0.4735900E-01,0.4743600E-01, + + 0.4752100E-01,0.4760900E-01,0.4769300E-01,0.4777400E-01, + + 0.4785200E-01,0.4791600E-01,0.4796000E-01,0.4799100E-01, + + 0.4801600E-01,0.4803900E-01,0.4806000E-01,0.4807950E-01, + + 0.2430000E-01,0.3474000E-01,0.3942000E-01,0.4135500E-01, + + 0.4203000E-01,0.4227300E-01,0.4235850E-01,0.4239270E-01, + + 0.4241520E-01,0.4243590E-01,0.4245660E-01,0.4247910E-01, + + 0.4250970E-01,0.4255920E-01,0.4262310E-01,0.4269240E-01, + + 0.4276890E-01,0.4284810E-01,0.4292370E-01,0.4299660E-01, + + 0.4306680E-01,0.4312440E-01,0.4316400E-01,0.4319190E-01, + + 0.4321440E-01,0.4323510E-01,0.4325400E-01,0.4327155E-01, + + 0.2255000E-01,0.3225500E-01,0.3659500E-01,0.3838900E-01, + + 0.3901500E-01,0.3924050E-01,0.3931985E-01,0.3935155E-01, + + 0.3937240E-01,0.3939160E-01,0.3941080E-01,0.3943165E-01, + + 0.3946005E-01,0.3950600E-01,0.3956530E-01,0.3962960E-01, + + 0.3970055E-01,0.3977400E-01,0.3984415E-01,0.3991180E-01, + + 0.3997695E-01,0.4003040E-01,0.4006715E-01,0.4009305E-01, + + 0.4011390E-01,0.4013310E-01,0.4015065E-01,0.4016695E-01, + + 0.2130000E-01,0.3044500E-01,0.3455500E-01,0.3625450E-01, + + 0.3684700E-01,0.3706050E-01,0.3713575E-01,0.3716575E-01, + + 0.3718550E-01,0.3720370E-01,0.3722190E-01,0.3724165E-01, + + 0.3726850E-01,0.3731195E-01,0.3736805E-01,0.3742890E-01, + + 0.3749605E-01,0.3756555E-01,0.3763190E-01,0.3769590E-01, + + 0.3775750E-01,0.3780805E-01,0.3784280E-01,0.3786725E-01, + + 0.3788700E-01,0.3790520E-01,0.3792180E-01,0.3793720E-01, + + 0.2025000E-01,0.2895000E-01,0.3285000E-01,0.3446250E-01, + + 0.3502500E-01,0.3522750E-01,0.3529875E-01,0.3532725E-01, + + 0.3534600E-01,0.3536325E-01,0.3538050E-01,0.3539925E-01, + + 0.3542475E-01,0.3546600E-01,0.3551925E-01,0.3557700E-01, + + 0.3564075E-01,0.3570675E-01,0.3576975E-01,0.3583050E-01, + + 0.3588900E-01,0.3593700E-01,0.3597000E-01,0.3599325E-01, + + 0.3601200E-01,0.3602925E-01,0.3604500E-01,0.3605960E-01, + + 0.1920000E-01,0.2745500E-01,0.3114500E-01,0.3267050E-01, + + 0.3320300E-01,0.3339470E-01,0.3346215E-01,0.3348915E-01, + + 0.3350690E-01,0.3352320E-01,0.3353950E-01,0.3355725E-01, + + 0.3358140E-01,0.3362045E-01,0.3367085E-01,0.3372550E-01, + + 0.3378585E-01,0.3384835E-01,0.3390800E-01,0.3396550E-01, + + 0.3402090E-01,0.3406635E-01,0.3409760E-01,0.3411965E-01, + + 0.3413740E-01,0.3415370E-01,0.3416860E-01,0.3418245E-01 / + + DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0,0.D0 / + DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0,1.D-9/ + DATA LAHG / 2000.0D2,4.0D5,1.0D6,4.0D6,1.0D7 / + DATA Lm /3.70D5,3.85D5,4.0D5,4.17D5,4.17D5,4.35D5,5.0D5,5.56D5, + + 5.99D5,6.33D5, 6.67D5,7.04D5 / + DATA LONG / + + 280.0, + + 300.0, + + 320.0, + + 340.0, + + 360.0, + + 380.0, + + 400.0, + + 450.0, + + 500.0, + + 550.0, + + 600.0, + + 650.0 / + DATA PI / 3.141592654D0 / + PARAMETER (RT=6348.0D5) + +*********************************************************************** +* +* SCATTERING PARAMETERS FOR RAYLEIGH: +* +* MEAN FREE PATH FOR SCATTERING RAYLEIGH (g/cm^2) + PARAMETER (XR=2.970D3) +*********************************************************************** + + + T_Ray = 1.0 + T_Mie = 1.0 + T_Oz = 1.0 + + H = -RT + SQRT(RT**2 + (height/COS(theta))**2 + + + (2.0D0*RT*height)) + + ROW = AINT(((H+1.)/1.0E5)) + +*********************************************************************** +* +* LARGE ZENITH ANGLE FACTOR (AIR MASS): + + + m= (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)- + + (RT*COS(theta))) +********************************************************************** + + +* RAYLEIGH SCATTERING + + RHOTOT = 0.0 + do 11 i=1,5 + RHOF(i) = 0 + 11 continue + + DO 91 I=2,5 + IF ( H .LT. LAHG(I) ) THEN + RHOTOT = RHOTOT + + + BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) - + + EXP(-H/CATM(I-1))) + GOTO 92 + ELSE + RHOTOT = RHOTOT + + + BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) - + + EXP(-LAHG(I)/CATM(I-1))) + ENDIF + 91 CONTINUE + + 92 RHO_FI = m*RHOTOT + + T_Ray = EXP(-(RHO_FI/XR)*(400.D0/wavelength)**4) + + +*************************************************************************** +* * +* MIE ABSORPTION: WE USE FOR HEIGHTS LOWER THAN 10 Km AN EXACT FORMULA* +* FOR THE AEROSOL DENSITY, AND WE USE A TABLE FOR HIGHTS HIGHER THAN * +* 10 Km. WE TAKE THE TABLE FROM THE ALTERMAN ARTICLE. * +*************************************************************************** + + +*************************************************************************** +* * +* MIE SCATTERING PARAMETERS * + Hm = 1.2D5 +*************************************************************************** + + IF (ROW.GT.27) THEN + ROW=28 + ENDIF + + + CON_MI = 2 + + 1001 IF (ABS(LONG(CON_MI)-wavelength).LT.0.01) THEN + TOT_AE =AE_ABI(ROW,CON_MI) + ELSEIF (LONG(CON_MI).GT.wavelength) THEN + A = (AE_ABI(ROW,CON_MI)-AE_ABI(ROW,CON_MI-1))/ + + (LONG(CON_MI) - LONG(CON_MI-1)) + B = AE_ABI(ROW,CON_MI) - (A*LONG(CON_MI)) + TOT_AE = A*wavelength + B + ELSE + CON_MI = CON_MI +1 + GOTO 1001 + ENDIF + T_Mie = EXP(-(m*TOT_AE)) + + + +*********************************************************************** +* * +* OZONE ABSORPTION: We used the Elterman table. * +* * +*********************************************************************** + IF (ROW.GT.47) THEN + ROW = 47 + ENDIF + + CON_OZ = 2 + 2001 IF (LONG(CON_OZ).GT.wavelength) THEN + A = (OZ_ABI(ROW,CON_OZ)-OZ_ABI(ROW,CON_OZ-1))/ + + (LONG(CON_OZ) - LONG(CON_OZ-1)) + B = OZ_ABI(ROW,CON_OZ) - (A*LONG(CON_OZ)) + TOT_OZ = (A* wavelength)+B + ELSEIF (ABS(LONG(CON_OZ)-wavelength).LT.0.01) THEN + TOT_OZ = OZ_ABI(ROW,CON_OZ) + ELSE + CON_OZ = CON_OZ + 1 + GOTO 2001 + ENDIF + + T_Oz = EXP(-(m*TOT_OZ)) + +************************************************************************ +* * +* TOTAL TRANSMISSION OF THE ATMOSPHERE: (RAYLEIGH + MIE + OZONE) * +************************************************************************ + + tr_atmos = T_Ray*T_Mie*T_Oz + + RETURN + + END + +c @endcode + +c EOF Index: trunk/MagicSoft/Simulation/Detector/Reflector/config.mk =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/config.mk (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/config.mk (revision 292) @@ -0,0 +1,61 @@ +################################################################## +# +# config.mk +# +# @file config.mk +# @title small configuration file for Makefile +# @author J C Gonz\'alez +# @email gonzalez@mppmu.mpg.de +# @date Fri Mar 12 11:51:11 MET 1999 +# +#_______________________________________________________________ +# +# Created: Fri Mar 12 11:51:11 MET 1999 +# Author: Jose Carlos Gonzalez +# Purpose: Makefile for the compilation of the reflector program +# Notes: +# +#--------------------------------------------------------------- +# $RCSfile: config.mk,v $ +# $Revision: 1.1.1.1 $ +# $Author: harald $ +# $Date: 1999-10-29 07:00:33 $ +################################################################## +# @maintitle + +# @code + +# program + +PROGRAM = camera + +# compilers & tools + +CC = cc +CXX = cxx +F77 = f77 + +DOCUM = ../sus/sus +RATE = ../../comstat/comstat -g -p + +# includes + +INCLUDE = ../include-GENERAL +INCLUDE_COR = ../include-CORSIKA +INCLUDE_MC = ../include-MC +INCLUDE_REFL = . + +OPTIM = -O +DEBUG = -g + +# libraries + +RANLIBDIR = ../lib +CERNDIR = /CERN + +# system + +SYSTEM = osf + +# @endcode +##EOF Index: trunk/MagicSoft/Simulation/Detector/Reflector/param =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/param (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/param (revision 292) @@ -0,0 +1,37 @@ +reflector 0.3 -*- sh -*- +# @T \newpage +# @section Source code of {\tt param} +# @code +# +# Sample parameters file +# +verbose_level 0 +# +#fixed_target 0. 0. +#range_events 1 2946 +#skip 1 2946 +#range_events 381 382 +#max_events 2 +#energy_cuts 800. 1000. +# +ct_file magic.def +# +#output_file /hd80/MC/CT1/ana/g_magic_run16.rfl +output_file ./p.dummy +# +#atm_model ATM_CORSIKA +atm_model ATM_NOATMOSPHERE +# +data_paths 1 +./ +# +data_from_stdin +data_to_stdout +# +#output_file ./pb.rfl +#parallel_beam 0. 0. 0. 0. 2000. 2000. 300 300 1.e7 +#pm_parallel_beam unnamed.ppm 2. 10000. +# +end_file +# +# @endcode Index: trunk/MagicSoft/Simulation/Detector/Reflector/readparam.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/readparam.cxx (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/readparam.cxx (revision 292) @@ -0,0 +1,853 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= readparam +//= +//= @file readparam.cxx +//= @desc Reading of parameters file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Program for reflector simulation +//= Notes: See files README for details +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: readparam.cxx,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-10-29 07:00:33 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |readparam.cxx|. + +/*!@" + + This section describes briefly the source code for the file + |readparam.cxx|. All the defines it uses are located in the files + |readparam.h| and |atm.h|. In the first one we can see the + definitions of the commands available for the parameters file. We + will describe these commands later on. Please, note that this + program is still in development, and more commands will be added + soon. + + @"*/ + +//!@subsection Includes and Global variables definition. + +/*!@" + + All the defines are located in the file |readparam.h|. + + @"*/ + +//!@{ + +#include "readparam.h" + +//!@} + +//!@subsection Definition of global variables. + +/*!@" + + Here we define the global variables where the values from the + parameters file are stored. + + @"*/ + +//!@{ + +//@: list of paths to get data from +static char **Paths_list; + +//@: number of paths +static int Num_of_paths = 1; + +//@: output filename +static char Output_filename[PATH_MAX_LENGTH]; + +//@: name of the CT def. file +static char CT_filename[PATH_MAX_LENGTH]; + +//@: verbosity +static VerboseLevel verbose = VERBOSE_DEFAULT; + +//@: fixed target zenith angle (theta) +static float fixed_Theta; + +//@: fixed target zenith angle (theta) +static float fixed_Phi; + +//@: lower limits for the energy of the showers to be used +static float low_Ecut = 0.0; + +//@: lower limits for the energy of the showers to be used +static float high_Ecut = 100000.0; + +//@: are we using fixed target? +static int is_Fixed_Target = FALSE; + +//@: maximum number of events to read +static int max_Events = 50000; + +//@: range of events to read +static int first_Event=0, last_Event=1000000; + +//@: random seeds +static long int Seeds[2] = {3141592L, 2718182L}; + +//@: flag: blocking? +static int Block=0; + +//@: flag: random pointing? +static int Random_Pointing = FALSE; + +//@: flag: random pointing maximum distance +static float Random_Pointing_MaxDist; + +//@: flag: do we read from STDIN? +static int Data_From_STDIN = FALSE; + +//@: flag: do we write to STDOUT? +static int Data_To_STDOUT = FALSE; + +//!@} + +/*!@" + + Here we define the global variables where the parameters of a + parallel beam of light will be stored. + + @"*/ + +//!@{ + +static int pb_ParallelBeam = FALSE; +static int pb_ParallelBeamPM = FALSE; +static float pb_Theta, pb_Phi, pb_X, pb_Y; +static float pb_Height; +static float pb_LengthX, pb_LengthY, pb_NX, pb_NY; +static float pb_Scale; +static char pb_Filename[40]; + +//!@} + +//!@subsection The function |readparam()|. + +//!--------------------------------------------------------------------- +// @name creadparam +// +// @desc read parameters from the stdin / parameters file +// +// @var *filename Name of the parameters file (NULL->STDIN) +// +// @date Mon Sep 14 13:27:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +void +readparam(char * filename) +{ + char sign[] = GLUE_postp( PROGRAM, VERSION ); // initialize sign + char line[LINE_MAX_LENGTH]; // line to get from the stdin + char token[ITEM_MAX_LENGTH]; // a single token + char atm[ATM_NAME_MAX_LENGTH]; // name of the atm. model from stdin + int i, j; // dummy counters + ifstream ifile; + + // use cin or ifile (reading from STDIN or from parameters file? + if ( filename != NULL ) + ifile.open( filename ); + + // get signature + if ( filename != NULL ) + ifile.getline(line, LINE_MAX_LENGTH); + else + cin.getline(line, LINE_MAX_LENGTH); + line[strlen(SIGNATURE)] = '\0'; + strcpy(line, sign); + if (strcmp(sign, SIGNATURE) != 0) { + cerr << "ERROR: Signature of parameters file is not correct\n"; + cerr << '"' << sign << '"' << '\n'; + cerr << "should be: " << SIGNATURE << '\n'; + exit(1); + } + + // loop till the "end" directive is reached + int is_end = FALSE; + while (! is_end) { + + // get line from file or stdin + if ( filename != NULL ) + ifile.getline(line, LINE_MAX_LENGTH); + else + cin.getline(line, LINE_MAX_LENGTH); + + // skip comments (start with '#') + if (line[0] == '#') + continue; + + // show user comments (start with '>') + if (line[0] == '>') { + cout << line << endl << flush; + continue; + } + + // look for each item at the beginning of the line + for (i=0; i<=end_file; i++) + if (strstr(line, ITEM_NAMES[i]) == line) + break; + + // if it is not a valid line, just ignore it + if (i == end_file+1) { + cerr << "Skipping unknown token in [" << line << "]\n"; + continue; + } + + // case block for each directive + switch ( i ) { + + case data_paths: // beginning of the list of valid paths + + // get number of paths to read + sscanf(line, "%s %d", token, &Num_of_paths); + + if ( verbose == TRUE ) + cout << Num_of_paths << " path(s) will be used.\n"; + + // allocate memory for paths list + Paths_list = (char**)malloc(Num_of_paths * sizeof(char*)); + for (i=0; i no blocking) +// +// @return Block value +// +// @date Wed Nov 25 13:21:00 MET 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +int +get_block(void) +{ + return ( Block ); +} +//!@} + + +//!--------------------------------------------------------------------- +// @name get_data_from_stdin +// +// @desc get whether we will read the data from the STDIN +// +// @return TRUE: we get data from STDIN; FALSE: oth. +// +// @date Wed Nov 25 13:21:00 MET 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +int +get_data_from_stdin(void) +{ + return ( Data_From_STDIN ); +} +//!@} + + +//!--------------------------------------------------------------------- +// @name get_data_to_stdout +// +// @desc get whether we will write the output data to the STDOUT +// +// @return TRUE: we put data into the STDOUT; FALSE: oth. +// +// @date Wed Nov 25 13:21:00 MET 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +int +get_data_to_stdout(void) +{ + return ( Data_To_STDOUT ); +} +//!@} + + +//!--------------------------------------------------------------------- +// @name get_random_pointing +// +// @desc get whether random pointing is selected +// +// @return TRUE: we use a random pointing +// +// @date Wed Nov 25 13:21:00 MET 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +int +get_random_pointing(float *maxdist) +{ + *maxdist = Random_Pointing_MaxDist; + return ( Random_Pointing ); +} +//!@} + + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ +// +// $Log: not supported by cvs2svn $ +// Revision 1.13 1999/10/05 11:06:38 gonzalez +// Sep. 1999 +// +// Revision 1.12 1999/03/24 16:33:00 gonzalez +// REFLECTOR 1.1: Release +// +// Revision 1.11 1999/01/21 16:03:37 gonzalez +// Only small modifications +// +// Revision 1.10 1999/01/19 18:07:18 gonzalez +// Bugs in STDIN-STDOUT version corrected. +// +// Revision 1.9 1999/01/14 17:35:46 gonzalez +// Both reading from STDIN (data_from_stdin) and +// writing to STDOUT (data_to_STDOUT) working. +// +// Revision 1.8 1999/01/13 12:41:11 gonzalez +// THIS IS A WORKING (last?) VERSION +// +// Revision 1.7 1998/12/15 10:47:24 gonzalez +// RELEASE-1.0 +// +//!@} + +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Reflector/readparam.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/readparam.h (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/readparam.h (revision 292) @@ -0,0 +1,230 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= readparam +//= +//= @file readparam.h +//= @desc Header file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Program for reflector simulation +//= Notes: See files README for details +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: readparam.h,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-10-29 07:00:33 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |readparam.h|. + +/*!@" + + In this section you can find the source code for the file + |readparam.h|. This file is mainly needed by |readparam.cxx|. Here + is located the definition of the commands you can use in the + parameters file. In this file, the first line must be |reflector| + |version|, where |version| is the appropiate version program which + can read the commands written in that parameters file. You cannot + mix parameters files and executables with different versions. The + current version is |reflector 0.3|. Each line starting with a \# + will be treated as a comment. + + The commands now available for the parameters file will be explained + in the documentation. + + A parameters file looks like this: + + |reflector 0.2| + + |verbose_level 2| + + |fixed_target 0. 0.| + + |max_events 1000| + + |data_paths 1| + + |/hd80/MC/CT1/data/low-za/p/run9| + + |output_file reflector.rfl| + + |ct_file ct1.def| + + |atm_model ATM_CORSIKA| + + |end_file| + + @"*/ + +//!@{ +#ifndef _readparam_ +#define _readparam_ + +#ifndef _this_ +#define _this_ readparam +#endif +//!@} + +//!@subsection Include files. + +//!@{ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "jcmacros.h" +#include "jcdebug.h" + +#include "reflector-v.h" + +#include "atm.h" + +//!@} + +//!@subsection Macro-definitions, and constants. + +//!@{ + +// now we define the list ITEM_LIST of possible items in +// the parameters file. note that they are functions of +// T(x). we will change T(x) to be the real item or the +// string of this item when needed + +#define ITEM_LIST /* LIST OF ITEMS IN THE PARAMETERS FILE */ \ +T(data_paths), /* begin of the list of source data paths */ \ +T(output_file), /* output file */ \ +T(ct_file), /* file with the characteristics of the CT */ \ +T(atm_model), /* changes the atmospheric model to be used */ \ +T(verbose_level), /* defines verbose level of the output */ \ +T(fixed_target), /* defines the position towards the CT is pointing */ \ +T(max_events), /* maximum number of event to read */ \ +T(range_events), /* analyze from event A to B */ \ +T(energy_cuts), /* lowest/highest energy allowed */ \ +T(parallel_beam), /* defines a parallel beam of light */ \ +T(pm_parallel_beam),/* defines a parallel beam of light, out of a pixmap */ \ +T(seeds), /* seeds for random number generation */ \ +T(block), /* size of the block of files, when 'blocking' */ \ +T(random_pointing), /* random CT pointing from each shower (hadrons) */ \ +T(data_from_stdin), /* read data from STDIN */ \ +T(data_to_stdout), /* read data from STDIN */ \ +T(end_file) /* end of the parameters file */ + +#define T(x) x // define T() as the name as it is + +enum ITEM_TYPE { + ITEM_LIST +}; + +#undef T + +#define T(x) #x // define T() as the string of x + +const char *const ITEM_NAMES[] = { + ITEM_LIST +}; + +#undef T + +#define LINE_MAX_LENGTH 400 +#define ITEM_MAX_LENGTH 40 +#define PATH_MAX_LENGTH 120 + +// Verbose Levels +enum VerboseLevel { + VERBOSE_QUIET, + VERBOSE_MINIMAL, + VERBOSE_NORMAL, + VERBOSE_MAXIMAL +}; + +#define VERBOSE_DEFAULT VERBOSE_QUIET + +// directory for intermediate storage of data from STDIN + +#define TMP_STDIN_DIR "/hd90/gonzalez/tmp" + +//!@} + +//!@subsection Prototypes of functions. + +//!@{ + +//++ +// prototypes +//-- + +void readparam(char * filename); +int get_num_of_paths(void); +char *get_path_name(int i); +char *get_output_filename(void); +char *get_ct_filename(void); +int get_verbose(void); +int get_fixed_target(float *theta, float *phi); +int get_max_events(void); +void get_range_events(int *a, int *b); +void get_energy_cuts(float *low, float *high); +void get_parallel_beam(float *pb_theta, float *pb_phi, + float *pb_x, float *pb_y, + float *pb_lengthx, float *pb_lengthy, + float *pb_nx, float *pb_ny, float *pb_height); +int is_parallel_beam(void); +char * get_parallel_beam_pm(float *pb_scale, float *pb_height); +int is_parallel_beam_pm(void); +long int get_seeds(int n); +int get_block(void); +int get_data_from_stdin(void); +int get_data_to_stdout(void); +int get_random_pointing(float *maxdist); + +//!@} + +//!@{ + +#endif // ! _readparam_ + +//!@} + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ + +/* + * $Log: not supported by cvs2svn $ + * Revision 1.13 1999/10/05 11:06:37 gonzalez + * Sep. 1999 + * + * Revision 1.12 1999/03/24 16:33:02 gonzalez + * REFLECTOR 1.1: Release + * + * Revision 1.11 1999/01/21 16:03:41 gonzalez + * Only small modifications + * + * Revision 1.10 1999/01/19 18:07:18 gonzalez + * Bugs in STDIN-STDOUT version corrected. + * + * Revision 1.9 1999/01/14 17:35:40 gonzalez + * Both reading from STDIN (data_from_stdin) and + * writing to STDOUT (data_to_STDOUT) working. + * + */ + +//!@} +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Reflector/reflector-v.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/reflector-v.h (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/reflector-v.h (revision 292) @@ -0,0 +1,12 @@ +#ifndef VERSION + +#define PROGRAM reflector +#define VERSION 0.3 + +#define GLUE_prep(x,y) #x" "#y +#define GLUE_postp(x,y) GLUE_prep(x,y) + +const char SIGNATURE[] = GLUE_postp( PROGRAM, VERSION ); + +#endif // ! VERSION + Index: trunk/MagicSoft/Simulation/Detector/Reflector/reflector.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/reflector.cxx (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/reflector.cxx (revision 292) @@ -0,0 +1,3567 @@ +//!///////////////////////////////////////////////////////////////////// +// +// reflector +// +// @file reflector.cxx +// @title Simulation of the Reflector +// @desc Program for the simulation of CT1 and MAGIC reflectors +// @author J C Gonzalez +// @email gonzalez@mppmu.mpg.de +// @date Thu May 7 16:24:22 1998 +// +//---------------------------------------------------------------------- +// +// Created: Thu May 7 16:24:22 1998 +// Author: Jose Carlos Gonzalez +// Purpose: Program for reflector simulation +// Notes: +// +//---------------------------------------------------------------------- +// +// $RCSfile: reflector.cxx,v $ +// $Revision: 1.1.1.1 $ +// $Author: harald $ +// $Date: 1999-10-29 07:00:33 $ +// +//////////////////////////////////////////////////////////////////////// +// @tableofcontents @coverpage + +//= Below you can find a sort of, incomplete for sure, documentation +//= concerning this program. Please, do not hesiate to contact the +//= author in case of problems. + +//=--------------------------------------------------------------------- +//!@section Details of the simulation of the detectors MAGIC and CT1. + +//=--------------------------------------------------------------------- +//!@subsection Overview of the simulation. + +/*!@" + + This short section gives an idea about how the simulation + developes itself. The details on each step will be explained + later. At the beginning we have the Generated output file. + + @itemize + + @- B@Atmospheric Absorption@: A detailed description of the + atmospheric absorption is provided. This description takes into + account the effects of the Rayleigh scattering, the aerosol (Mie) + absorption, and the effect of the Ozone. It works for the wavelength + range @$290 \div 600@$ nm. + + @- B@Reflection in the mirrors@: Each remaining photon is traced + till the mirrors, and then reflected till the camera plane. The + steps are the following: + + @enumerate + + @- Calculation of the mirror where the photon is + hitting. + + @- We apply the E@reflectivity@ effect, using a continuous random + number in the range [0:1] : the photon can be lost. + + @- If the photon is hitting in the E@black spot@ of the center, + the photon is also lost. + + @- We take the E@mathematicaly exact@ normal of the mirror, and + apply the E@axis-deviation@ effect: the normal changes a little + bit. + + @- We apply the E@imperfections@ in the mirrors (measurements are + needed). + + @- Finally, with the displaced normal, we calculate the position of + the reflected photon in the camera plane. Here, even with an exact, + mathematical ray-tracing the image is E@blurred@, due to the + intrinsic behaviour of the telescope optics. + + @endenumerate + + At this step we have a B@RAW-Monte-Carlo@ output file. + + @- B@Collection in the PMTs/IPCs@: The collection in the active + devices follows the next steps: + + @enumerate + + @- We apply a certain probability of transmision of the photon in + the E@plexiglas@, dependent on the incident angle of the photon in + the camera (measurements are needed). + + @- We apply the transmitance factor of the E@light guides@. + + \item Then we apply the E@Quantum efficiency@ of the PMT, + depending on the incident angle. + + @endenumerate + + @- B@Correction of the obtained signal@ In every pixel we have a + certain number of Cherenkov photons, what we take as a B@mean value@ + of the measured charge in that pixel. Together to that value, we + take into account the E@Night Sky Background@ (NSB): we have for + this a mean value measured, and with this mean value we calculate + the contribution of the sky for each pixel. + + At this level we will have a B@RAW-Data@ output file. + + @enditemize + + @"*/ + +//=--------------------------------------------------------------------- +//!@subsection Simulation of the detector. + +//=--------------------------------------------------------------------- +//!@subsubsection Coordinate systems. + +/*!@" + + The different angles used in the simulation of the detector are + shown in the next figures. The telescope (CT) is pointing to the + direction @$(\theta,\phi)@$, while the particle (background) is coming + from the direction @$(\theta_p,\phi_p)@$. This direction is calculated + using a pair of random generated angles @$\theta_s@$ and @$\phi_s@$, and + spherical trigonometry. + + @F + \begin{figure}[htbp] + \begin{center} + \includegraphics[height=0.45\textheight]{coordinates.eps} + \caption{Coordinate system used in the reflector simulation} + \label{fig:coord} + \end{center} + \end{figure} + @F + + @F + \stepcounter{figure} + \begin{figure}[htbp] + \begin{center} + \includegraphics[height=0.45\textheight]{reflector.eps} + \caption{Angles used in the reflector simulation} + \label{fig:reflector} + \end{center} + \end{figure} + @F + + The angles are defined in the next table. Using this definitions, we + can start to explain how the reflector program works. + + @F + \begin{table} + \begin{tabular}[hbt]{cl} + $\theta$ & Zenith angle of the position where the CT points to\\ + $\phi$ & Horizontal angle from the North of the position\\ + & where the CT points to\\ + $\theta_s$ & Angular distance between the axis of the\\ + & incoming shower and the direction where the CT is \\ + & pointing to.\\ + $\phi_s$ & Polar angle with axis equal to the axis of the CT.\\ + $\theta_p$ & Zenith angle of the incoming shower.\\ + $\phi_p$ & Horizontal angle (from the North) of the\\ + & direction of the incoming shower.\\ + $\theta_{CT}$ & Zenith angle of the dish ($= \theta$).\\ + $\phi_{CT}$ & Horizontal angle of the dish ($= \phi$).\\ + $\theta_m$ & Angle of the mirror axis relative to the CT axis.\\ + $\phi_m$ & Polar angle of the position of the mirror + (arbitrary origin).\\ + \end{tabular} + \caption{Definition of the different angles in the reflector program.} + \end{table} + @F + + At the CORSIKA level, and also at the reflector simulation level, we + calculate the angles @$\theta_p@$ and @$\phi_p@$, which give the + direction of the incoming shower. To calculate them we use the + following. + + @F + \begin{eqnarray} + \cos\theta_p &=& \cos\theta \cos\theta_s + + \sin\theta \sin\theta_s \cos\phi_s + \\ + \sin\phi_p &=& \frac{\sin\theta_s}{\sin\theta_p} \sin\phi_s + \end{eqnarray} + @F + + To finish with the definitions, we just write down some of + the parameters we will use in the reflexion algorithm, in + the next table. + + @F + \stepcounter{table} + \begin{table} + \begin{tabular}[hbt]{cl} + $\mathbf{x}=(x,y,z)$ & + Position of the Cherenkov photon on the ground\\ + & in the global system OXYZ.\\ + + $\mathbf{r}=(u,v,w)$ & + Vector of the Cherenkov photon's trayectory\\ + & in the global system OXYZ.\\ + + $\mathbf{x_{CT}}=(x_{CT},y_{CT},z_{CT})$ & + Position of the Cherenkov photon on the ground\\ + & in the system OXYZ$_{CT}$ of the CT.\\ + + $\mathbf{r_{CT}}=(u_{CT},v_{CT},w_{CT})$ & + Vector of the Cherenkov photon's trayectory\\ + & in the system OXYZ$_{CT}$ of the CT.\\ + + $\mathbf{x_{m}}=(x_{m},y_{m},z_{m})$ & + Position of the Cherenkov photon on the ground\\ + & in the system OXYZ$_{m}$ of the mirror.\\ + + $\mathbf{r_{m}}=(u_{m},v_{m},w_{m})$ & + Vector of the Cherenkov photon's trayectory\\ + & in the system OXYZ$_{m}$ of the mirror.\\ + + $\mathbf{x_{cut}}=(x_{cut},y_{cut},z_{cut})$ & + Cut of the trayectory of the photon with the mirror,\\ + & in the system OXYZ$_{m}$ of the mirror.\\ + + $\mathbf{r_\perp}=(u_\perp,v_\perp,w_\perp)$ & + Vector perpendicular to the mirror in the point $\mathbf{x_{cut}}$\\ + & in the system OXYZ$_{m}$ of the mirror.\\ + + $\mathbf{r_r}=(u_r,v_r,w_r)$ & + Vector of the reflected photon, \\ + & in the system OXYZ$_{m}$ of the mirror.\\ + + $\mathbf{r_r^{CT}}=(u_r^{CT},v_r^{CT},w_r^{CT})$ & + Vector of the reflected photon, \\ + & in the system OXYZ$_{CT}$ of the CT.\\ + + $\mathbf{x_{cam}}=(x_{cam},y_{cam},z_{cam})$ & + Position of the photon in the camera plane,\\ + & in the system OXYZ$_{CT}$ of the CT.\\ + + \end{tabular} + \caption{Different vectors and points used in the reflector + program.} + \end{table} + @F + + @"*/ + +//=--------------------------------------------------------------------- +//!@subsubsection Algorithm of the reflexion. + +/*!@" + + We start calculating the matrices for changing of system of + coordinates, @$\Omega(\theta,\phi)@$ and + @$\Omega^{-1}(\theta,\phi)@$. They are defined to be: + + @F + \begin{equation} + \begin{split} + \Omega(\theta,\phi) &= R_y(-\theta) R_z(-\phi) \\ + \Omega^{-1}(\theta,\phi) &= R_z(\phi) R_y(\theta) + \end{split} + \end{equation} + @F + + where the matrices @$R_x@$, @$R_y@$ and @$R_z@$ are rotations of + angle @$\alpha@$ around each axis OX, OY and OZ. + + @F + \begin{eqnarray} + R_x(\alpha) = + \begin{bmatrix} + 1 & 0 & 0 \\ + 0 & \cos\alpha & -\sin\alpha \\ + 0 & \sin\alpha & \cos\alpha \\ + \end{bmatrix} + \label{eq:rotx} + \\ + R_y(\alpha) = + \begin{bmatrix} + \cos\alpha & 0 & \sin\alpha \\ + 0 & 1 & 0 \\ + -\sin\alpha & 0 & \cos\alpha \\ + \end{bmatrix} + \label{eq:roty} + \\ + R_z(\alpha) = + \begin{bmatrix} + \cos\alpha & -\sin\alpha & 0\\ + \sin\alpha & \cos\alpha & 0\\ + 0 & 0 & 1 \\ + \end{bmatrix} + \label{eq:rotz} + \end{eqnarray} + @F + + With this, the matrices @$\Omega(\theta,\phi)@$ and + @$\Omega^{-1}(\theta,\phi)@$ result to be: + + @F + \begin{eqnarray} + \Omega(\theta,\phi) = + \begin{bmatrix} + \cos\phi \cos\theta&\sin\phi \cos\theta&-\sin\theta \\ + -\sin\phi & \cos\phi & 0 \\ + \cos\phi \sin\theta&\sin\phi \sin\theta& \cos\theta \\ + \end{bmatrix} + \\ + \Omega^{-1}(\theta,\phi) = + \begin{bmatrix} + \cos\phi \cos\theta&-\sin\phi&\cos\phi \sin\theta\\ + \sin\phi \cos\theta& \cos\phi&\sin\phi \sin\theta\\ + -\sin\theta & 0 & \cos\theta \\ + \end{bmatrix} + \end{eqnarray} + @F + + and each change of system of coordinates uses these + matrices. In order to calculate the new coordinates X' of a + point X in the new system O', and viceversa, we apply the + following equations: + + @F + \begin{equation} + \begin{split} + X' &= \Omega(\theta,\phi) X\\ + X &= \Omega^{-1}(\theta,\phi) X'\\ + \end{split} + \end{equation} + @F + + The steps taken in the reflector simulation are the following: + + @enumerate + + @- We look for the mirror @$m@$ which is closer to the + trayectory of the photon (presumably, the photon will hit + it). + + @- We move to the system @$\mbox{OXYZ}_{CT}@$ of the CT. + + @- We move again, now to the system @$\mbox{OXYZ}_{m}@$ of the + mirror @$m@$. + + @- Then the point where the trayectory of the photon intersects with + the mirror is calculated. + + @- In this point, the normal (perpendicular) to the mirror is + obtained. + + @- The reflected trayectory is calculated with the normal in the + mirror and the original trayectory. + + @- We move back to the system @$\mbox{OXYZ}_{CT}@$, and calculate + there the point where the new trayectory of the photon hits the + camera plane. + + @endenumerate + + In between the previous points the effects of the atmospheric + absorption and those mentioned in the next section are included. At + the end, the position of the photon and the incident angle in the + camera plane is saved. + + @"*/ + +//=--------------------------------------------------------------------- +//!@subsection Effects to be included in the simulation of the reflector. + +/*!@" + + @itemize + + @- B@Incidence in the mirrors@: + For each single photon we trace the point where it hits a + given mirror. We can use either hexagonal mirrors or + circular mirrors with the same surface of the hexagonal ones. + + @- B@Reflectivity of the mirrors@: + The reflectivity of the mirrors is a function @$R(\lambda)@$ + of the wavelength @$\lambda@$ of the photon. + + @- B@Black spots@: + In the center of each mirror there is a ``bad'' region, + where no reflection can be done at all. The radii of these + regions vary, but they can be @$r \simeq 2 \mathrm{\,cm}@$. + + @- B@Imperfections@: + The surface of each mirror is locally not-perfect. We need + measurements to include this effect. + + @- B@Axis--deviation@: + The axis of each mirror can be deviated from the + ``teorical'' point in the center of the camera. (@$\sigma + \simeq 2.5 \mathrm{\,mm}@$). + + @- B@Blurring@: + Pure effect of optics: each mirror is ``seen'' by the camera + with a different angle. + + @enditemize + + @"*/ + +//=--------------------------------------------------------------------- +//!@subsection Effects to be included in the simulation of the camera. + +/*!@" + + @itemize + + @- B@Plexiglas@: + In front of the light guides there is a ``plexiglas'' layer, + with a given transmitance (dependent on the incidence angle). + + @- B@Light Guides@: + The looses in the transmitance in the light guides are + around @$15-20@$\%. + + @- B@Quantum Efficiency@: + The QE of each PMT depends on the impact point of the photon + in the window of the PMT and the angle of incidence in this + point. + + @- B@Excess noise factor@: + It produces a smearing of the signal with a sigma + @$\sqrt{2N}@$ instead of the typical @$\sqrt{N}@$. + + @- B@Pick-up noise?@: + + @- B@First dynode collection efficiency?@: + + @enditemize + + @"*/ + +//=--------------------------------------------------------------------- +// @T \newpage + +//=--------------------------------------------------------------------- +//!@section Source code of |reflector.cxx|. + +/*!@" + + In this section we show the (commented) code of the program + for the simulation of the reflector, in the current version. + + @"*/ + +//=--------------------------------------------------------------------- +//!@subsection Preliminary notes in the code. + +//!@{ +//================================================================= +// +// NOTE: Most of the equations used in this program can be found +// in Bronstein, Semendjajew, Musiol and Muehlig, "Taschenbuch +// der Mathematik", Verlag Harri Deutsch, Thun und Frankfurt +// am Main, 1997. Some others were "deduced". +// +// +// Matrices for rotation of coordinate systems +//--------------------------------------------- +// +// Omega(theta,phi) = Ry(theta) Rz(phi) +// Omega-1(theta,phi) = Rz(-phi) Ry(-theta) +// +// X` = Omega(theta,phi) X +// X = Omega-1(theta,phi) X` +// +// Omega(theta,phi) = +// / \ +// | cos(phi)cos(theta) sin(phi)cos(theta) -sin(theta) | +// | -sin(phi) cos(phi) 0 | +// | cos(phi)sin(theta) sin(phi)sin(theta) cos(theta) | +// \ / +// +// OmegaI(theta,phi) = +// / \ +// | cos(phi)cos(theta) -sin(phi) cos(phi)sin(theta) | +// | sin(phi)cos(theta) cos(phi) sin(phi)sin(theta) | +// | -sin(theta) 0 cos(theta) | +// \ / +// +//================================================================= +//!@} + +//=--------------------------------------------------------------------- +//!@subsection Includes and Global variables definition. + +/*!@" + + All the defines are located in the file |reflector.h|. + + @"*/ + +//!@{ +#include "reflector.h" +//!@} + +//=--------------------------------------------------------------------- +//!@subsubsection Definition of global variables. + +/*!@" + + We first define the matrices to change of coordinates system. The + definition of the matrices was already seen in the previous section. + All matrices are in the form @$m[i][j]@$ with i the running index + for rows, and j the one for columns. + + @"*/ + +//!@{ + +//@: matrices to change to the system where the optical axis is OZ +static float OmegaCT[3][3]; + +//@: matrices to change to the system where the optical axis is OZ (inverse) +static float OmegaICT[3][3]; + +//@: matrices to change the system of coordinates +static float Omega[3][3]; + +//@: matrices to change the system of coordinates (inverse) +static float OmegaI[3][3]; + +//!@} + +/*!@" + + Now we define some global variables with data about the telescope, + such as E@focal distance@, E@number of pixels/mirrors@, E@size of + the camera@, and so on. + + @"*/ + +/*!@" + + Depending on the telescope we are using (CT1 or MAGIC), the + information stored in the definition file is different. The + variable |ct_Type| has the value 0 when we use CT1, and 1 when we + use MAGIC. + + @"*/ + +//!@{ + +//@: Type of telescope (CT1: 0; MAGIC: 1) +static int ct_Type; + +//!@} + +/*!@" + + And this is the information about the whole telescope. + +@"*/ + +//!@{ +// parameters of the CT (from the CT definition file) + +//@: Focal distances [cm] +static float *ct_Focal; + +//@: Mean Focal distances [cm] +static float ct_Focal_mean; + +//@: STDev. Focal distances [cm] +static float ct_Focal_std; + +//@: Mean Point Spread function [cm] +static float ct_PSpread_mean; + +//@: STDev. Point Spread function [cm] +static float ct_PSpread_std; + +//@: STDev. Adjustmente deviation [cm] +static float ct_Adjustment_std; + +//@: Radius of the Black Spot in mirror [cm] +static float ct_BlackSpot_rad; + +//@: Radius of one mirror [cm] +static float ct_RMirror; + +//@: Camera width [cm] +static float ct_CameraWidth; + +//@: Pixel width [cm] +static float ct_PixelWidth; + +//@: Number of mirrors +static int ct_NMirrors = 0; + +//@: Number of pixels +static int ct_NPixels; + +//!@} + +/*!@" + + The following double-pointer is a 2-dimensional table with + information about each mirror in the dish. The routine + |read_ct_file()| will read this information from the file with the + name given by the user in the E@parameters file@, in the command + |ct_file|. The information stored in this file (and in this table) + depends on the type of telescope we are using. + + @"*/ + +//!@{ + +//@: Pointer to a table with the following info.: +static float **ct_data; + +/* + * TYPE=0 (CT1) + * i s rho theta x y z thetan phin xn yn zn + * + * i : number of the mirror + * s : arc length [cm] + * rho : polar rho of the position of the center of the mirror [cm] + * theta : polar angle of the position of the center of the mirror [cm] + * x : x coordinate of the center of the mirror [cm] + * y : y coordinate of the center of the mirror [cm] + * z : z coordinate of the center of the mirror [cm] + * thetan : polar theta angle of the direction where the mirror points to + * phin : polar phi angle of the direction where the mirror points to + * xn : xn coordinate of the normal vector in the center (normalized) + * yn : yn coordinate of the normal vector in the center (normalized) + * zn : zn coordinate of the normal vector in the center (normalized) + * + * TYPE=1 (MAGIC) + * i f sx sy x y z thetan phin + * + * i : number of the mirror + * f : focal distance of that mirror + * sx : curvilinear coordinate of mirror's center in X[cm] + * sy : curvilinear coordinate of mirror's center in X[cm] + * x : x coordinate of the center of the mirror [cm] + * y : y coordinate of the center of the mirror [cm] + * z : z coordinate of the center of the mirror [cm] + * thetan : polar theta angle of the direction where the mirror points to + * phin : polar phi angle of the direction where the mirror points to + * xn : xn coordinate of the normal vector in the center (normalized) + * yn : yn coordinate of the normal vector in the center (normalized) + * zn : zn coordinate of the normal vector in the center (normalized) + */ + +//!@} + +/*!@" + + Here we define the global variables where the parameters of a + parallel beam of light will be stored. + + @"*/ + +//!@{ + +static int pb_ParallelBeam = FALSE; +static int pb_ParallelBeamPM = FALSE; +static float pb_Theta, pb_Phi, pb_X, pb_Y; +static float pb_LengthX, pb_LengthY, pb_NX, pb_NY; +static float pb_Scale; +static float pb_Height; +static char pb_Filename[40]; + +//!@} + +/*!@" + + We can analyze the data stored in the directories provided via the + |data_paths| entry in the parameters file in one go, or in blocks of + size |Block|, if this variable is a positive integer. This can be + done through a shell script. The mechanism is the following: + + @itemize + + @- First, we read the parameters file, and wait for a file called + |__DOIT| (written by the script when a new block of |Block| files + is available in a loop. + + @- If this file exists, we remove it and create another file with + the name |__WORKING| (in order to let the script to know that we + are working), and process a block of |Block| files. After that, we + will rename the file |__WORKING| to |__READY|. + + @- When the script sees a new block of |Block| files, it will write + again a file called |__DOIT|, unless there exist a file called + |__WORKING| (which means that the |reflector| is still working in + the previous block). + + @- Each time a block is processed, is customary to the user to + delete or not that block of files (or to save it to tape, ...). + + @- At the end, a file called |__DONE| is written telling us the + number of blocks processed. + + @enditemize + + @"*/ + +//!@{ + +static int Block=0; +int WORKING=FALSE; +int nfilesBlock=0; +ofstream obfile; +long posdir; + +//!@} + +/*!@" + + We can read also the data from the standard input, and we could + write data into the standard output. + + @"*/ + +//!@{ + +//@: boolean variable to control whether we read from STDIN or not +static int Data_From_STDIN = FALSE; + +//@: boolean variable to control whether we write to STDOUT or not +static int Data_To_STDOUT = FALSE; + +//!@} + +/*!@" + + Next, we have two tables with data got from two files, + |reflectivity.dat| and |axisdev.dat|, which contain information + about the reflectivity of the mirrors, as a function of the + wavelength @$\lambda@$, and about the (simulated) deviation of the + mirrors' axes, with respect to the mathematical exact axes. + + @"*/ + +//!@{ + +// table of reflectivity for each mirror + +//@: table with datapoints (wavelength,reflec.) +static float **Reflectivity; + +//@: number of datapoints +static int nReflectivity; + +//@: table with deviations of the mirrors' normals +static float **AxisDeviation; + +//!@} + +/*!@" + + We still define a table into where normal random numbers will be + stored by the routine |rnormal(double *r, int n)|. + + @"*/ + +//!@{ + +//@: table of normal random numbers +static double NormalRandomNumbers[500]; + +//!@} + +/*!@" + + This is a flag to change the E@verbosity@ of the output + + @"*/ + +//!@{ + +//@: flag to change the verbosity +int verbose = VERBOSE_DEFAULT; + +//!@} + +/*!@" + + This option makes the telescope to point to a random position + relative to the shower axis, with a maximum angular separation of + |RANDOM_POINTING_MAX_SEPARATION| (defined in |reflector.h|. + + @"*/ + +//!@{ + +//@: random pointing for the CT? +int Random_Pointing = FALSE; + +//@: maximum random pointing distance +float Random_Pointing_MaxDist; // [radians] + +//!@} + +//=--------------------------------------------------------------------- +//!@subsection Main program. + +/*!@" + + In the main program we basically follow the items expressed in the + first section. Please, have a look into previous sections. I tried + to make a code already well commented. Most things should be clear. + Nevertheless, I'll try to comment again everything in following + releases. + + @"*/ + +//!@{ + +//++++++++++++++++++++++++++++++++++++++++ +// MAIN PROGRAM +//---------------------------------------- + +int +main(int argc, char **argv) +{ + + DIR *directory; //@< directory of data + struct dirent *de; //@< directory entry + + char pathname[256]; //@< name of directory + char cername[256]; //@< name of Cherenkov file + char staname[256]; //@< name of Statistics file + char outname[256]; //@< name of the output file + char parname[256]; //@< name of parameters file + + ifstream cerfile; //@< stream for Cherenkov file + ifstream pmfile; //@< stream for pixmap file + ofstream outputfile; //@< stream for output file + ofstream pbfile; //@< stream for parallel beam file + + COREventHeader evth; //@< Event Header class (from CORSIKA) + CORParticle photon; //@< Particle (photon) class (from CORSIKA) + CORStatfile stat; //@< Statistics file class (from CORSIKA) + + MCEventHeader mcevth; //@< Event Header class + MCCphoton cphoton; //@< Particle (photon) class + + float thetaCT, phiCT, xiCT; //@< parameters of a given shower + float coreD, coreX, coreY; //@< core position and distance + float u, v, w; //@< cosenos directores + + float r[3]; //@< trayectory + float x[3]; //@< position of the photon in the ground + + float rCT[3]; //@< trayectory in the system of the CT + float xCT[3]; //@< position of the photon in the ground (CT) + float rm[3]; //@< trayectory in the system of a mirror + float xmm[3]; //@< intermediate values + float xm[3]; //@< position of photon in ground + float xcut[3]; //@< location of the cut sphere-trayectory + float xcutCT[3]; //@< location of the cut sphere-trayectory (CT) + + float rnor[3], rnorm; //@< normal in that point + + float rrefl[3]; //@< reflected vector, from the mirror + float rreflCT[3]; //@< reflected vector, from the CT + float xcam[3]; //@< where the photon hits the camera plane + + float calpha; //@< cos(alpha=angle incident/normal) + float phi; //@< angle between photon and camera plane + + float a, b, c, t; //@< intermediate variables + + float d; //@< minimum distance trayectory-mirror center + + float wl; //@< wavelength of the Cphoton + float reflec; //@< reflectivity for a Cphoton + + int ncer; //@< number of the shower; + int np; //@< number of path + int i, k; //@< simple counters + int i_mirror=-1; //@< number of a given mirror + + int nCphotons; //@< number of a Cphotons written on output file + + float TR; //@< atmospheric transmitance + int nbeforeTR, nafterTR; //@< number of Cph. before and after transmission + int num_cer_files; //@< number of cer files; + int max_num_cer_files; //@< maximum number of cer files to read + int f_ev, l_ev; //@< first and last events of a range + + float lE, uE; //@< lower and upper edges of energy range allowed + + float distmirr, distmirr2; //@< distances used in MAGIC reflexion routine + float sx, sy; + float t1, t2; + float pbx, pby, pbt, pbp; + + char line[200]; + + int ch, errflg = 0; //@< used by getopt + + /*!@' + + We start in the main program. First we (could) make some + presentation, and follows the reading of the parameters file (now + from the STDIN), the reading of the CT parameters file, and the + creation of the output file, where the processed data will be + stored. + + */ + + //++ + // START + //-- + + // make unbuffered output + + cout.setf ( ios::stdio ); + + // parse command line options (see reflector.h) + + parname[0] = '\0'; + + optarg = NULL; + while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) ) + switch (ch) { + case 'f': + strcpy(parname, optarg); + break; + case 'h': + usage(); + break; + default : + errflg++; + } + + // show help if error + + if ( errflg>0 ) + usage(); + + // make some sort of presentation + // present(); + + // read parameters from the stdin + + if ( strlen(parname) < 1 ) + readparam(NULL); + else + readparam(parname); + + // blocking-mode? + + Block = get_block(); + + // reading data from STDIN? + + Data_From_STDIN = get_data_from_stdin(); + + // writing data to STDOUT? + + Data_To_STDOUT = get_data_to_stdout(); + + // make it verbosely? + + verbose = get_verbose(); + + /*!@' + + @#### Parallel beam of light. + + In order to test the behaviour of the program, we can feed it with + an artificial data file. This data file will be generated by the + program itself and it will consist in a grid of photons, forming + all together a parallel beam of light with certain + characteristics. These characteristics will be given by the user + in the |parameters file|, in the following form. + + parallel\_beam | | + + or + + parallel\_beam\_pm | | + + In the first form, the first two values give the arrival direction + of the beam (zenith angle and azimuth) in degrees, and the last + six values specify the geometry of the beam, in cm, and the height + of production, also in cm. In the second form, we use a pixmap + file in PPM ASCII format, and scale it to create (for 0 degrees of + zenith angle) a parallel beam of light. + + */ + + //++ + // PARALLEL BEAM OF LIGHT + //-- + + // are we going to use a parallel beam of light? + + pb_ParallelBeam = is_parallel_beam(); + pb_ParallelBeamPM = is_parallel_beam_pm(); + + if ( pb_ParallelBeam == TRUE ) { + + if ( pb_ParallelBeamPM == FALSE ) { + + // get parameters of the parallel beam of light + + get_parallel_beam(&pb_Theta, &pb_Phi, + &pb_X, &pb_Y, + &pb_LengthX, &pb_LengthY, + &pb_NX, &pb_NY, &pb_Height); + + // generate artificial file(s) + + // a. cerfile + + sprintf( cername, "%s/cer999999", get_path_name(0) ); + + pbfile.open( cername, ios::out ); + + evth.fill ( 1., 1., 1000., 1., 1.e6, 0., 0., 0., pb_Theta, pb_Phi, + pb_X, pb_Y); + + evth.write( pbfile ); + + for ( pbx = (pb_X - pb_LengthX/2.0); + pbx <= (pb_X + pb_LengthX/2.0); + pbx += (pb_LengthX/pb_NX) ) + + for ( pby = (pb_Y - pb_LengthY/2.0); + pby <= (pb_Y + pb_LengthY/2.0); + pby += (pb_LengthY/pb_NY) ) { + + pbt = -atan2(sqrt(SQR(pbx)+SQR(pby)),pb_Height); + pbp = atan2(pby, pbx); + photon.fill(450.0, + pbx, pby, + -sin(pbt) * cos(pbp), + -sin(pbt) * sin(pbp), + 0., pb_Height); + + photon.write( pbfile ); + + } + + photon.fill( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ); + + photon.write( pbfile ); + + pbfile.close(); + + // b. statfile + + sprintf( staname, "%s/sta999999", get_path_name(0) ); + + pbfile.open( staname ); + + stat.write( pbfile ); + + pbfile.close(); + + } else { + + // get parameters of the parallel beam of light from PIXMAP + + strcpy(pb_Filename, get_parallel_beam_pm(&pb_Scale,&pb_Height)); + + pmfile.open( pb_Filename, ios::out ); + pmfile.getline( line, 300 ); + + do { + pmfile.getline( line, 300 ); + } while (line[0] == '#'); + + sscanf( line, "%f %f", &pb_NX, &pb_NY ); + pb_LengthX = pb_NX * pb_Scale; + pb_LengthY = pb_NY * pb_Scale; + + // generate artificial file(s) + + // a. cerfile + + sprintf( cername, "%s/cer999999", get_path_name(0) ); + + pbfile.open( cername, ios::out ); + + evth.fill ( 1., 1., 1000., 1., 1.e6, 0., 0., 0., 0., 0., 0., 0.); + + evth.write( pbfile ); + + for ( pbx = (pb_X - pb_LengthX/2.0); + pbx <= (pb_X + pb_LengthX/2.0); + pbx += (pb_LengthX/pb_NX) ) + + for ( pby = (pb_Y - pb_LengthY/2.0); + pby <= (pb_Y + pb_LengthY/2.0); + pby += (pb_LengthY/pb_NY) ) { + + pmfile >> i; + + if ( i < 1 ) + continue; + + pbt = -atan2(sqrt(SQR(pbx)+SQR(pby)),pb_Height); + pbp = atan2(pby, pbx); + photon.fill(450.0, + pbx, pby, + -sin(pbt) * cos(pbp), + -sin(pbt) * sin(pbp), + 0., pb_Height); + + photon.write( pbfile ); + + } + + photon.fill( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ); + + photon.write( pbfile ); + + pbfile.close(); + + // b. statfile + + sprintf( staname, "%s/sta999999", get_path_name(0) ); + + pbfile.open( staname ); + + stat.write( pbfile ); + + pbfile.close(); + + pmfile.close(); + + } + + } // YES parallel beam + + // read parameters from the stdin + + read_ct_file(); + + // set all random numbers seeds + + setall( get_seeds(0), get_seeds(1) ); + + // get name of the output file + + strcpy( outname, get_output_filename() ); + + // open((re)create) output file + + if ( verbose >= VERBOSE_MINIMAL ) + log( SIGNATURE, "Openning/creating output file %s\n", outname ); + + outputfile.open( outname ); + + if ( outputfile.bad() ) + error( SIGNATURE, "Cannot open output file: %s\n", outname ); + + // save signature to the output file + + if ( Data_To_STDOUT ) + cout.write( SIGNATURE, sizeof(SIGNATURE) ); + else + outputfile.write( SIGNATURE, sizeof(SIGNATURE) ); + + // generate a sort of log information + + if ( verbose >= VERBOSE_MINIMAL ) { + log( SIGNATURE, "Atmospheric model: %s\n", get_atm_mod() ); + log( SIGNATURE, "Number of paths: %d\n", get_num_of_paths() ); + log( SIGNATURE, "Output filename: %s\n", outname ); + log( SIGNATURE, "CT def. filename: %s\n", get_ct_filename() ); + } + + // get energy range + + get_energy_cuts( &lE, &uE ); + + /*!@' + + After reading the parameters file (from STDIN), and the CT + definition file, we start the analysis. Basically, there are three + loops. The outermost loop is a loop in the E@data directories@ we + introduced in the parameters file; the middle loop follows each + |cer*| file in the directory; and the innermost loop gets each + Cherenkov from these files. + + */ + + // for each path of data files + + for (np=0; np= VERBOSE_MINIMAL ) + log( SIGNATURE, "Openning directory %s\n", pathname ); + + directory = opendir(pathname); + + if ( directory == NULL ) + error( SIGNATURE, + "Cannot open directory %s\n", pathname ); + + // write start-of-run mark + + if ( Data_To_STDOUT ) + cout.write( FLAG_START_OF_RUN, SIZE_OF_FLAGS ); + else + outputfile.write( FLAG_START_OF_RUN, SIZE_OF_FLAGS ); + + // trace the number of events (cer* files) + + num_cer_files = 0; + + // get maximum number, first and last event + + max_num_cer_files = get_max_events(); + + get_range_events(&f_ev, &l_ev); + + // if reading from disk, save position of handler in directory + + posdir = telldir( directory ); + + if ( Data_From_STDIN ) { + + if ( verbose >= VERBOSE_MINIMAL ) + log(SIGNATURE, " * * * READING DATA FROM STDIN * * *\n"); + + get_stdin_files(NULL, lE, uE, TRUE); + + } + + //--- LOOP ON FILES --- + // for each directory entry (files), + // and(or) while still data has to be read + // this condition is: + // while ( ( there are files AND + // maximum number is not reached ) + // OR + // ( reading from stdin AND + // maximum number is not reached ) ) + + while ( ( + ((de = readdir( directory )) != NULL) && + (num_cer_files < max_num_cer_files) + ) + || + ( + ( Data_From_STDIN ) && + (num_cer_files < max_num_cer_files) + ) ) { + + // if Block > 0, then we wait till a file __DOIT is present + + if ( Block > 0 ) { + + if ( ! WORKING ) { + + while ( ! WORKING ) { + // sleep 10 minutes + sleep( 60 ); + + // look for the file __DOIT + obfile.open("__DOIT"); + if ( ! obfile.bad() ) { + WORKING = TRUE; + obfile.close(); + } + } + + obfile.open("__WORKING"); + remove("__DOIT"); + + rewinddir( directory ); + seekdir( directory, posdir ); + + nfilesBlock = 0; + + continue; + + } else { + + if ( nfilesBlock == Block ) { + + posdir = telldir( directory ); + obfile.close(); + obfile.open("__READY"); + obfile << nfilesBlock << " files/block, " + << ncer << " showers.\n"; + obfile.close(); + remove("__WORKING"); + WORKING = FALSE; + + continue; + } + + } + + } // if Block + + if ( Data_From_STDIN ) { + + if ( verbose >= VERBOSE_MINIMAL ) + log(SIGNATURE, "Try to get data from STDIN\n"); + + if ( ! get_stdin_files(num_cer_files) ) { + if ( verbose >= VERBOSE_MINIMAL ) + log(SIGNATURE, "get_stdin_files returned FALSE.\n"); + if (cin.eof()) { + num_cer_files = max_num_cer_files + 1; + continue; + } + } + + ++num_cer_files; + ncer = num_cer_files; + + } else { + + // skip removed files + + if ( de->d_ino == 0) + continue; + + // keep only cer* files + + if ( strstr(de->d_name, "cer") != de->d_name ) + continue; + + if ( Block > 0 ) { + ++nfilesBlock; + obfile << nfilesBlock; + } + + // get cer* file number (number of the shower) + + // increments the pointer in 3 to calculate the number + + ncer = atoi(de->d_name + 3); + + if ( (ncer < f_ev) || (ncer > l_ev) ) + continue; + + // it is a real cer* file, and we want to read it + + ++num_cer_files; + + if ( verbose >= VERBOSE_NORMAL ) + log( SIGNATURE, "cerfile: %s\n", de->d_name ); + + } // data from disk + + // Here we initialize the counters of Cherenkov photons at + // different levels of the analysis. + + nCphotons = nbeforeTR = nafterTR = 0; + + // full names of the input files cer* and sta* + + sprintf(cername, "%s/cer%06d", pathname, ncer); + sprintf(staname, "%s/sta%06d", pathname, ncer); + + // try to open cer* file + + if ( verbose >= VERBOSE_MAXIMAL ) + log( SIGNATURE, "Opening %s\n", cername ); + + cerfile.open( cername ); + + if ( cerfile.bad() ) + error( SIGNATURE, + "Cannot open input file: %s\n", cername ); + + /*!@' + + Each shower has associated three files: E@ Cherenkov-photons + file @ (|cer*|), E@ Particles file @ (|dat*|) and E@ + Statistics file @ (|sta*|). First we read the data from this + last file, which is explained below, and then we read the E@ + Event-Header @ block from the |cer*| file. + + */ + + // try to open the statistics file + + if ( verbose >= VERBOSE_MAXIMAL ) + log( SIGNATURE, "Opening %s\n", staname ); + + stat.openfile ( staname ); + + // reading data from this sta* file + + if ( verbose >= VERBOSE_MAXIMAL ) + log( SIGNATURE, "Reading %s\n", staname ); + + stat.read(); + + // read the Event Header from the cer* file + + evth.read( cerfile ); + + if ( verbose >= VERBOSE_MAXIMAL ) { + log( SIGNATURE, "EVENT HEADER:\n"); + evth.print(); + } + + // check energy + + if ( (evth.get_energy() < lE) || + (evth.get_energy() > uE) ) { + + // close files + cerfile.close(); + stat.closefile(); + + // skip it + if ( verbose >= VERBOSE_MAXIMAL ) + log( SIGNATURE, "Bad energy.\n"); + + continue; + } + + // write start-of-event mark + + if ( Data_To_STDOUT ) + cout.write( FLAG_START_OF_EVENT, SIZE_OF_FLAGS ); + else + outputfile.write( FLAG_START_OF_EVENT, SIZE_OF_FLAGS ); + + // save information from the cer* to the mc-evth + + mcevth.transport( &evth ); + + // save time interval in the mc-evth + + mcevth.put_times ( stat.get_tfirst(), stat.get_tlast() ); + + // get direction where the CT is pointing to + // (or, better, from when the shower is coming from) + + // first, put the values from the shower + thetaCT = evth.get_theta(); + phiCT = evth.get_phi(); + + // then, may be we have to change them + + // do we want random pointing (around shower axis) ? + if ( get_random_pointing( &Random_Pointing_MaxDist ) == TRUE ) { + + // we do, then get a random position + xiCT = get_new_ct_pointing(thetaCT, phiCT, + Random_Pointing_MaxDist, + &thetaCT, &phiCT); + + xiCT = xiCT; // thi is to avoid "was set but not used" messages + + // save deviations in mc-evth + mcevth.put_deviations (thetaCT - evth.get_theta(), + phiCT - evth.get_phi()); + + if ( verbose >= VERBOSE_MAXIMAL ) + log( SIGNATURE, "NEW POINTING: (%f,%f) -> (%f,%f)\n", + evth.get_theta(), evth.get_phi(), thetaCT, phiCT); + + } else { + + mcevth.put_deviations (0.0, 0.0); + + if ( get_fixed_target( &thetaCT, &phiCT ) == FALSE ) { + thetaCT = evth.get_theta(); + phiCT = evth.get_phi(); + } + + } + + // save event-header to the output file + + if ( Data_To_STDOUT ) + cout.write ( (char *)&mcevth, mcevth.mysize() ); + else + mcevth.write( outputfile ); + + if ( verbose >= VERBOSE_MAXIMAL ) + mcevth.print(); + + // get core distance and position + + coreD = evth.get_core(&coreX, &coreY); + coreD=coreD; // to avoid "was set but not used" messages + + // calculate matrices for changing of coordinate system from-to CT + + makeOmega( thetaCT, phiCT ); + makeOmegaI( thetaCT, phiCT ); + + memcpy( OmegaCT, Omega, 9*sizeof(float) ); + memcpy( OmegaICT, OmegaI, 9*sizeof(float) ); + + // read each and every single photon (particle) from cer* file + + while ( ! (cerfile.eof() || cerfile.bad()) ) { + + // read photon data from the file + + photon.read( cerfile ); + + wl = photon.get_wl(); + + if ( wl < 1.0 ) + break; + + // get director cosines in z-axis + + w = r[2] = photon.get_w(); + + /*!@' + + @#### Atmospheric attenuation. + + The routine |atm()| calls the proper attenuation routine, + and returns the transmission coefficient for that + whavelength, height and zenith angle. + + */ + + //++ + // FILTER No. 1: ATMOSPHERE + //-- + + ++nbeforeTR; + + // get transmitance of the atmosphere + + TR = atm( photon.get_wl(), photon.get_h(), acos(w) ); + + // passes the atmosphere? + + // if random > TR then goes to the TOP of the loop again + + if ( RandomNumber > TR ) + continue; + + ++nafterTR; + + /*!@' + + @#### Reflectivity of the mirrors. + + We make a 3rd. order interpolation using Lagrange + polynomials, in order to calculate the reflectivity of the + mirror for that wavelength. Further developments will + include also a incidence-angle dependence (which is not very + important). + + */ + + //++ + // FILTER No. 2: REFLECTIVITY R(lambda) + //-- + + // find data point to be used in Lagrange interpolation (-> k) + + FindLagrange(Reflectivity,k,wl); + + // if random > reflectivity then goes to the TOP of the loop again + + reflec = Lagrange(Reflectivity,k,wl); + + if ( RandomNumber > reflec ) + continue; + + /*!@' + + @#### Reflexion in the mirrors. + + The reflexion is made now is a somehow dirty way: if the + value of the variable |ct\_Type| is 0, the telescope is CT1, + and if the value is 1, the telescope is MAGIC. Due to the + completly different geometries of the dishes and the + mirrors, now we support two fully separated and different + algorithms. + + */ + + //++ + // REFLEXION + //-- + + // get director cosines x,y and save original trayectory + + u = r[0] = photon.get_u(); + v = r[1] = photon.get_v(); + u=u, v=v; // to avoid "was set but not used" messages + + // get position in the ground + + x[0] = photon.get_x() - coreX; + x[1] = photon.get_y() - coreY; + x[2] = 0.0; + + if ( verbose >= VERBOSE_MAXIMAL ) { + cout << '@' << '1' + << ' ' << ncer + << ' ' << t + << ' ' << x[0] + << ' ' << x[1] + << ' ' << x[2] + << ' ' << r[0] + << ' ' << r[1] + << ' ' << r[2] + << ' ' << thetaCT + << ' ' << phiCT + << endl; + } + + // change to the system of the CT + + applyMxV( OmegaCT, x, xCT ); + applyMxV( OmegaCT, r, rCT ); + + // before moving to the system of the mirror, for MAGIC, + // first we look whether the photon hits a mirror or not + + // calculate the intersection of the trayectory of the photon + // with the GLOBAL DISH !!! + // we reproduce the calculation of the coefficients of the + // second order polynomial in z (=xCT[2]), made with + // Mathematica + + /* + * In[1]:= parab:=z-(x^2+y^2)/(4F) + * par1=parab /. {x->x0+u/w(z-z0),y->y0+v/w(z-z0)} + * + * Out[1]= + * u (z - z0) 2 v (z - z0) 2 + * (x0 + ----------) + (y0 + ----------) + * w w + * z - --------------------------------------- + * 4 F + * + * In[2]:= CoefficientList[ExpandAll[par1*4F*w^2],z] + * + * Out[2]= + * 2 2 2 2 + * {-(w x0 ) - w y0 + 2 u w x0 z0 + + * + * 2 2 2 2 + * 2 v w y0 z0 - u z0 - v z0 , + * + * 2 2 + * 4 F w - 2 u w x0 - 2 v w y0 + 2 u z0 + + * + * 2 2 2 + * 2 v z0, -u - v } + */ + + // the z coordinate is calculated + + a = - SQR(rCT[0]) - SQR(rCT[1]); + b = 4.0*ct_Focal_mean*SQR(rCT[2]) + - 2.0*rCT[0]*rCT[2]*xCT[0] - 2.0*rCT[1]*rCT[2]*xCT[1] + + 2.0*SQR(rCT[0])*xCT[2] + 2.0*SQR(rCT[1])*xCT[2]; + c = 2*rCT[0]*rCT[2]*x[0]*x[2] + 2*rCT[1]*rCT[2]*x[1]*x[2] + - SQR(rCT[2])*SQR(x[0]) - SQR(rCT[2])*SQR(x[1]) + - SQR(rCT[0])*SQR(x[2]) - SQR(rCT[1])*SQR(x[2]); + + if ( fabs(a) < 1.e-6 ) { + + // only one value + + xcut[2] = -c / b; + + } else { + + d = sqrt( b*b - 4.0*a*c ); + + // two possible values for z + + t1 = (-b+d) / (2.0*a); + t2 = (-b-d) / (2.0*a); + + // z must be the minimum of t1 and t2 + + xcut[2] = (t1 < t2) ? t1 : t2; + + } + + // xcut[] is NOW the cut between the GLOBAL dish of MAGIC and + // the trayectory of the photon + + xcut[0] = xCT[0] + rCT[0]/rCT[2]*(xcut[2]-xCT[2]); + xcut[1] = xCT[1] + rCT[1]/rCT[2]*(xcut[2]-xCT[2]); + + if ( ct_Type == 1 ) { + + // ++ >>>>>>>>>> MAGIC <<<<<<<<<< + + // convert to Curvilinear distance over the parabolic dish + + sx = Lin2Curv( xcut[0] ); + sy = Lin2Curv( xcut[1] ); + + // is it outside the dish? + + if ((fabs(sx) > 850.0) || (fabs(sy) > 850.0)) { + //cout << "CONDICION 1 !" << endl << flush; + //cout << '1'; + continue; + } + + // -- >>>>>>>>>> MAGIC <<<<<<<<<< + + } // endif (ct_Type == 1) + + // calculate the mirror to be used + + distmirr = 1000000.; + + for (i=0; i=ct_RMirror; ++i) { + distmirr2 = sqrt(SQR(ct_data[i][CT_X] - xcut[0]) + + SQR(ct_data[i][CT_Y] - xcut[1]) + + SQR(ct_data[i][CT_Z] - xcut[2])); + if (distmirr2 < distmirr) { + i_mirror = i; + distmirr = distmirr2; + } + } + + // the mirror to use is i_mirror (calculated several lines above) + // check whether the photon is outside the nearest (this) mirror + + if ( ct_Type == 0 ) { + + // ++ >>>>>>>>>> CT1 <<<<<<<<<< + if (distmirr > ct_RMirror) { + //cout << "CONDICION 2 !" << endl << flush; + //cout << '2'; + continue; + } + // -- >>>>>>>>>> CT1 <<<<<<<<<< + + } else { + + // ++ >>>>>>>>>> MAGIC <<<<<<<<<< + if ((fabs(ct_data[i_mirror][CT_SX] - sx) > ct_RMirror) || + (fabs(ct_data[i_mirror][CT_SY] - sy) > ct_RMirror)) { + //cout << "CONDICION 2 !" << endl << flush; + //cout << '2'; + continue; + } + // -- >>>>>>>>>> MAGIC <<<<<<<<<< + + } + + // calculate matrices for the mirror + + makeOmega (-ct_data[i_mirror][CT_THETAN], + ct_data[i_mirror][CT_PHIN]); + makeOmegaI(-ct_data[i_mirror][CT_THETAN], + ct_data[i_mirror][CT_PHIN]); + + // change to the system of the mirror + + // first translation... + + xmm[0] = xCT[0] - ct_data[i_mirror][CT_X]; + xmm[1] = xCT[1] - ct_data[i_mirror][CT_Y]; + xmm[2] = xCT[2] - ct_data[i_mirror][CT_Z]; + + // ...then rotation + + applyMxV( Omega, xmm, xm ); + applyMxV( Omega, rCT, rm ); + + // the vector rCT should be normalized, and + // so the vector rm remains normalized as well, but, anyhow... + + rnorm = NORM( rm ); + rm[0] /= rnorm; + rm[1] /= rnorm; + rm[2] /= rnorm; + + // calculate the intersection of the trayectory of the photon + // with the mirror + // we reproduce the calculation of the coefficients of the + // second order polynomial in z (=xm[2]), made with + // Mathematica + + /* + * In[1]:= esfera:=x^2+y^2+(z-R)^2-R^2; + * recta:={x->x0+u/w(z-z0),y->y0+v/w(z-z0)} + * + * In[2]:= esfera + * + * 2 2 2 2 + * Out[2]= -R + x + y + (-R + z) + * + * In[3]:= recta + * + * u (z - z0) v (z - z0) + * Out[3]= {x -> x0 + ----------, y -> y0 + ----------} + * w w + * + * In[4]:= esf=esfera /. recta + * + * 2 2 u (z - z0) 2 v (z - z0) 2 + * Out[4]= -R + (-R + z) + (x0 + ----------) + (y0 + ----------) + * w w + * + * In[5]:= coefs=CoefficientList[ExpandAll[esf],z] + * + * 2 2 2 2 + * 2 2 2 u x0 z0 2 v y0 z0 u z0 v z0 + * Out[5]= {x0 + y0 - --------- - --------- + ------ + ------, + * w w 2 2 + * w w + * + * 2 2 2 2 + * 2 u x0 2 v y0 2 u z0 2 v z0 u v + * > -2 R + ------ + ------ - ------- - -------, 1 + -- + --} + * w w 2 2 2 2 + * w w w w + * In[6]:= Simplify[ExpandAll[coefs*w^2]] + * + * 2 2 2 2 2 2 + * Out[6]= {w (x0 + y0 ) - 2 w (u x0 + v y0) z0 + (u + v ) z0 , + * + * 2 2 2 2 2 + * > -2 (R w - u w x0 + u z0 + v (-(w y0) + v z0)), u + v + w } + * + */ + + // the z coordinate is calculated, using the coefficients + // shown above + + a = SQR(rm[0]) + SQR(rm[1]) + SQR(rm[2]); + if ( ct_Type == 0 ) { + // CT1 + b = -2*(2.*ct_Focal_mean*SQR(rm[2]) + - rm[0]*rm[2]*xm[0] + + SQR(rm[0])*xm[2] + + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2])); + } else { + // MAGIC + b = -2*(2.*ct_data[i_mirror][CT_FOCAL]*SQR(rm[2]) + - rm[0]*rm[2]*xm[0] + + SQR(rm[0])*xm[2] + + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2])); + } + c = (SQR(rm[2])*(SQR(xm[0]) + SQR(xm[1])) + - 2*rm[2]*(rm[0]*xm[0] + rm[1]*xm[1])*xm[2] + + (SQR(rm[0]) + SQR(rm[1]))*SQR(xm[2])); + + d = sqrt( b*b - 4.0*a*c ); + + // two possible values for z + + t1 = (-b+d) / (2.0*a); + t2 = (-b-d) / (2.0*a); + + // z must be the minimum of t1 and t2 + + xcut[2] = (t1 < t2) ? t1 : t2; + xcut[0] = xm[0] + rm[0]/rm[2]*(xcut[2]-xm[2]); + xcut[1] = xm[1] + rm[1]/rm[2]*(xcut[2]-xm[2]); + + //++ + // BLACK SPOTS: If the photon hits the black spot, it's lost + //-- + + if ( sqrt(SQR(xcut[0]) + SQR(xcut[1])) < ct_BlackSpot_rad ) { + //cout << "CONDITION 3!\n" << flush; + //cout << '3'; + continue; + } + + // if still we have the photon, continue with the reflexion + + // calculate normal vector in this point + // (and normalize, with the sign changed) + + rnor[0] = 2.0*xcut[0]; + rnor[1] = 2.0*xcut[1]; + rnor[2] = 2.0*(xcut[2] - 2.0*ct_Focal[i_mirror]); + + rnorm = -NORM( rnor ); + rnor[0] /= rnorm; + rnor[1] /= rnorm; + rnor[2] /= rnorm; + + // now, both "normal" vector and original trayectory are + // normalized + // just project the original vector in the normal, and + // take it as the "mean" position of the original and + // the "reflected" vector + // from this, we can calculate the "reflected" vector + // calpha = cos(angle(rnor,rm)) + + calpha = fabs(rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]); + + // finally!!! we have the reflected trayectory of the photon + + rrefl[0] = 2.0*rnor[0]*calpha - rm[0]; + rrefl[1] = 2.0*rnor[1]*calpha - rm[1]; + rrefl[2] = 2.0*rnor[2]*calpha - rm[2]; + + rnorm = NORM( rrefl ); + rrefl[0] /= rnorm; + rrefl[1] /= rnorm; + rrefl[2] /= rnorm; + + // let's go back to the coordinate system of the CT + + // first rotation... + applyMxV( OmegaI, xcut, xcutCT); + applyMxV( OmegaI, rrefl, rreflCT); + + // ...then translation + xcutCT[0] += ct_data[i_mirror][CT_X]; + xcutCT[1] += ct_data[i_mirror][CT_Y]; + xcutCT[2] += ct_data[i_mirror][CT_Z]; + + // calculate intersection of this trayectory and the camera plane + // in the system of the CT, this plane is z = ct_Focal + + t = (ct_Focal_mean - xcutCT[2]) / rreflCT[2]; + + xcam[0] = xcutCT[0] + rreflCT[0]*t; + xcam[1] = xcutCT[1] + rreflCT[1]*t; + xcam[2] = xcutCT[2] + rreflCT[2]*t; + + //++ + // AXIS DEVIATION: We introduce it here just as a first order + // correction, by modifying the position of the reflected photon. + //-- + + xcam[0] += AxisDeviation[0][i_mirror]; + xcam[1] += AxisDeviation[1][i_mirror]; + + //++ + // SMEARING: We apply the point spread function for the mirrors + //-- + + // get two N(0;1) random numbers + + rnormal( NormalRandomNumbers, 2 ); + + // modify the Cphoton position in the camera + + xcam[0] += NormalRandomNumbers[0] * ct_PSpread_mean; + xcam[1] += NormalRandomNumbers[1] * ct_PSpread_mean; + + // check whether the photon goes out of the camera + + if ( (SQR(xcam[0])+SQR(xcam[1])) > SQR(ct_CameraWidth) ) { + continue; + } + + //++ + // ANGLE OF INCIDENCE + //-- + + // calculate angle of incidence between tray. and camera plane + // the camera plane is + // 0 y + 0 y + z - ct_Focal = 0 => (A,B,C,D) = (0,0,1,-ct_Focal) + // from Table 3.20 "Tasch. der Math." + + phi = asin(rreflCT[2]); + + //++ + // TIMING + //-- + + // calculate the new time of the photon (in the camera) + + // initial time since first interaction + + t = photon.get_t(); + + // substract path from the mirror till the ground, 'cos + // the photon actually hit the mirror!! + + t = t - sqrt( SQR(xm[0] - xcut[0]) + + SQR(xm[1] - xcut[1]) + + SQR(xm[2] - xcut[2]) ) / Speed_of_Light_air_cmns; + + // add path from the mirror till the camera + + t = t + sqrt( SQR(xcutCT[0] - xcam[0]) + + SQR(xcutCT[1] - xcam[1]) + + SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns; + + // save data in the photon variable + + cphoton.fill( photon.get_id(), + xcam[0], xcam[1], + photon.get_u(), photon.get_v(), + photon.get_h(), t, phi); + + // show it + + if ( verbose >= VERBOSE_MAXIMAL ) { + cout << '@' << '2' + << ' ' << ncer + << ' ' << xCT[0] + << ' ' << xCT[1] + << ' ' << xCT[2] + << ' ' << rCT[0] + << ' ' << rCT[1] + << ' ' << rCT[2] + << ' ' << xcut[0] + << ' ' << xcut[1] + << ' ' << xcut[2] + << endl; + cout << '@' << '3' + << ' ' << ncer + << ' ' << sx + << ' ' << sy + << ' ' << i_mirror + << ' ' << ct_data[i_mirror][CT_SX] + << ' ' << ct_data[i_mirror][CT_SY] + << ' ' << ct_data[i_mirror][CT_SX] - sx + << ' ' << ct_data[i_mirror][CT_SY] - sy + << endl; + if ( pb_ParallelBeam == FALSE ) { + cout << '@' << '4' + << ' ' << ncer + << ' ' << xcam[0] + << ' ' << xcam[1] + << ' ' << xcam[2] + << ' ' << phi + << ' ' << photon.get_t()-stat.get_tfirst() + << ' ' << t-stat.get_tfirst() + << endl << flush; + } else { + cout << '@' << '4' + << ' ' << ncer + << ' ' << xcam[0] + << ' ' << xcam[1] + << ' ' << xcam[2] + << ' ' << phi + << endl << flush; + } + } + + // write the photon data to the output file + + if ( Data_To_STDOUT ) + cout.write ( (char *)&cphoton, cphoton.mysize() ); + else + cphoton.write( outputfile ); + + // increase number of Cphotons written + + ++nCphotons; + + } // while still there are photons left + + if ( verbose >= VERBOSE_NORMAL ) + cout << nCphotons << ' ' + << nbeforeTR << ' ' + << nafterTR << ' ' + << endl << flush; + + // write end-of-event mark + + if ( Data_To_STDOUT ) + cout.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS ); + else + outputfile.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS ); + + // close files + + cerfile.close(); + stat.closefile(); + + // show how many photons were written + + if ( verbose >= VERBOSE_MAXIMAL ) + log( SIGNATURE, "%d C-photons written.\n", nCphotons ); + + // if we are reading data from STDIN, delete the last used files + + if ( Data_From_STDIN ) { + + unlink( cername ); + unlink( staname ); + + } + + } // while there are still data left + + // write end-of-run mark + + if ( Data_To_STDOUT ) + cout.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS ); + else + outputfile.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS ); + + // close directory + + closedir( directory ); + + } // for each data directory + + // write end-of-file mark + + if ( Data_To_STDOUT ) + cout.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS ); + else + outputfile.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS ); + + // close output file + + if ( verbose >= VERBOSE_MINIMAL ) + log( SIGNATURE, "Closing output file %s\n", outname ); + + if ( ! Data_To_STDOUT ) + outputfile.close(); + + // clean everything + + if ( verbose >= VERBOSE_MINIMAL ) + log( SIGNATURE, "Cleanning . . .\n" ); + + clean(); + + // program finished + + if ( verbose >= VERBOSE_MINIMAL ) + log( SIGNATURE, "Done.\n"); + + return ( 0 ); +} +//!@} + +// @T \newpage + +//=--------------------------------------------------------------------- +// @subsection Functions definition + +//!--------------------------------------------------------------------- +// @name present +// +// @desc Make some presentation +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +void +present(void) +{ + cout << "##################################################\n" + << SIGNATURE << '\n' << '\n' + << "Simulation of the reflector\n" + << "J C Gonzalez, May 1998\n" + << "##################################################\n\n"; +} +//!@} + + +//!--------------------------------------------------------------------- +// @name usage +// +// @desc show help +// +// @date Tue Dec 15 16:23:30 MET 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +void +usage(void) +{ + present(); + cout << "\nusage ::\n\n" + << "\t reflector " + << " [ -@ paramfile ] " + << " [ -h ] " + << "\n\n or \n\n" + << "\t reflector < paramfile" + << "\n\n"; + exit(0); +} +//!@} + + +//!--------------------------------------------------------------------- +// @name clean +// +// @desc basically, frees memory +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +void +clean(void) +{ + int i; + + // delete focals table + + delete [] ct_Focal; + ct_Focal = NULL; + + // delete reflectivity table + + for (i=0; i<2; i++) { + delete [] Reflectivity[i]; + } + + delete [] Reflectivity; + Reflectivity = NULL; + + // delete mirrors' data table + + for (i=0; i= VERBOSE_MINIMAL ) + log( "read_ct_file", "start.\n" ); + + ifstream ctin ( get_ct_filename() ); + + if ( ctin.bad() ) + error( "read_ct_file", + "Cannot open CT def. file: %s\n", get_ct_filename() ); + + // loop till the "end" directive is reached + + while (!ctin.eof()) { + + // get line from stdin + + ctin.getline(line, LINE_MAX_LENGTH); + + // look for each item at the beginning of the line + + for (i=0; i<=define_mirrors; i++) + if (strstr(line, CT_ITEM_NAMES[i]) == line) + break; + + // if it is not a valid line, just ignore it + + if (i == define_mirrors+1) + continue; + + // case block for each directive + + switch ( i ) { + + case type: // (0:CT1 ¦ 1:MAGIC) + + // get focal distance + + sscanf(line, "%s %d", token, &ct_Type); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", ": %s\n", + ((ct_Type==0) ? "CT1" : "MAGIC") ); + + break; + + case focal_distance: // [cm] + + // get focal distance + + sscanf(line, "%s %f", token, &ct_Focal_mean); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", ": %f cm\n", ct_Focal_mean ); + + break; + + case focal_std: // s(focal distance) [cm] + + // get focal distance + + sscanf(line, "%s %f", token, &ct_Focal_std); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std ); + + break; + + case point_spread: // [cm] + + // get point spread + + sscanf(line, "%s %f", token, &ct_PSpread_mean); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", ": %f cm\n", ct_PSpread_mean ); + + break; + + case point_std: // s(point spread) [cm] + + // get point spread + + sscanf(line, "%s %f", token, &ct_PSpread_std); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std ); + + break; + + case adjustment_dev: // s(adjustment_dev) [cm] + + // get point spread + + sscanf(line, "%s %f", token, &ct_Adjustment_std); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std ); + + break; + + case black_spot: // radius of the black spot in the center [cm] + + // get black spot radius + + sscanf(line, "%s %f", token, &ct_BlackSpot_rad); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Radius of the black spots: %f cm\n", + ct_BlackSpot_rad); + + break; + + case r_mirror: // radius of the mirrors [cm] + + // get radius of mirror + + sscanf(line, "%s %f", token, &ct_RMirror); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror ); + + break; + + case n_mirrors: // number of mirrors + + // get the name of the output_file from the line + + sscanf(line, "%s %d", token, &ct_NMirrors); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors ); + + break; + + case camera_width: // camera width [cm] + + // get the name of the ct_file from the line + + sscanf(line, "%s %f", token, &ct_CameraWidth); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth ); + + break; + + case n_pixels: // number of pixels + + // get the name of the output_file from the line + + sscanf(line, "%s %d", token, &ct_NPixels); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels ); + + break; + + case pixel_width: // pixel width [cm] + + // get the name of the ct_file from the line + + sscanf(line, "%s %f", token, &ct_PixelWidth); + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth ); + + break; + + case define_mirrors: // read table with the parameters of the mirrors + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Table of mirrors data:\n" ); + + // check whether the number of mirrors was already set + + if ( ct_NMirrors == 0 ) + error( "read_ct_file", "NMirrors was not set.\n" ); + + // allocate memory for paths list + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_ct_file", "Allocating memory for ct_data\n" ); + + ct_data = new float*[ct_NMirrors]; + + for (i=0; i= VERBOSE_MINIMAL ) + log( "read_ct_file", "Reading mirrors data...\n" ); + + // read ASCII data + + for (i=0; i= VERBOSE_MINIMAL ) + log( "read_ct_file", "done.\n" ); + return; +} +//!@} + + +//!--------------------------------------------------------------------- +// @name dist_r_P +// +// @desc distance straight line r - point P +// +// @var a coord. X of a fixed point of the straight line +// @var b coord. Y of a fixed point of the straight line +// @var c coord. Z of a fixed point of the straight line +// @var l component X of a vector of the straight line +// @var m component Y of a vector of the straight line +// @var n component Z of a vector of the straight line +// @var x coord. X of the point P +// @var y coord. Y of the point P +// @var z coord. Z of the point P +// +// @return Returns the distance d(r,P) between the line and the point P +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +float +dist_r_P(float a, float b, float c, + float l, float m, float n, + float x, float y, float z) +{ + return ( + sqrt((SQR((a-x)*m-(b-y)*l) + + SQR((b-y)*n-(c-z)*m) + + SQR((c-z)*l-(a-x)*n))/ + (SQR(l)+SQR(m)+SQR(n)) + ) + ); +} +//!@} + + +//!--------------------------------------------------------------------- +// @name read_reflectivity +// +// @desc read reflectivity data for the mirrors +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +void +read_reflectivity(void) +{ + ifstream reflfile; + ofstream reflfile_write; + int i, errors=0; + char line[LINE_MAX_LENGTH]; + + // try to open the file + reflfile.open( REFLECTIVITY_FILE ); + + // if it is wrong or does not exist: + + while ( reflfile.bad() ) { + + ++errors; + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_reflectivity", "Cannot open reflectivity file: %s\n", + REFLECTIVITY_FILE ); + + if ( errors > 5 ) + error( "read_reflectivity", "Exiting."); + + // try to re-generate the file + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_reflectivity", "Generating file: %s\n", + REFLECTIVITY_FILE ); + + // try to open the file + reflfile_write.open( REFLECTIVITY_FILE ); + + // write short comment in the beginning + + reflfile_write << "# reflectivity file" << endl + << "# reflectivity for each mirror in the frame" << endl + << "# J C Gonzalez, 1998" << endl + << "#" << endl; + + // write data + + // write number of data-points + + reflfile_write << "# number of datapoints" << endl + << 35 << endl; + + // write data-points + + reflfile_write << "# datapoints" << endl; + + for ( i=0; i<35; ++i ) + reflfile_write << 270.0+i*10.0 << ' ' << .9 << endl; + + // close it and try to re-open + + reflfile_write.close(); + reflfile.close(); + reflfile.open( REFLECTIVITY_FILE ); + + } + + if ( verbose >= VERBOSE_MINIMAL ) { + if ( errors > 0 ) + log( "read_reflectivity", "Reading reflectivity file after %d fail%c\n", + errors, ((errors>1)?'s':' ') ); + else + log( "read_reflectivity", "Reading reflectivity file . . .\n" ); + } + + // we set this to -1, in order to read in the first place the + // number of data-points for the Reflectivity curve + i = -1; + + // scan the whole file + + while ( ! reflfile.eof() ) { + + // get line from the file + + reflfile.getline(line, LINE_MAX_LENGTH); + + // skip if comment + + if ( *line == '#' ) + continue; + + // get the value + + if ( i < 0 ) { // if it is the first number we read + + // read the number of datapoints + + sscanf(line, "%d", &nReflectivity); + + // allocate memory for Reflectivity table + + Reflectivity = new float * [2]; + Reflectivity[0] = new float[nReflectivity]; + Reflectivity[1] = new float[nReflectivity]; + + } else { // if not, then it is a datapoint + + // get the datapoint (wavelength, reflectivity) + + sscanf(line, "%f %f", &Reflectivity[0][i], &Reflectivity[1][i]); + + } + + ++i; + + } + + // close it + + reflfile.close(); + +} +//!@} + + +//!--------------------------------------------------------------------- +// @name read_axisdev +// +// @desc read axis deviations data for the mirrors +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +void +read_axisdev(void) +{ + ifstream axisdevfile; + ofstream axisdevfile_write; + int i, errors=0; + char line[LINE_MAX_LENGTH]; + + // try to open the file + axisdevfile.open( AXISDEVIATION_FILE ); + + // if it is wrong or does not exist: + + while ( axisdevfile.bad() ) { + + ++errors; + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_axisdev", "Cannot open axis deviation file: %s\n", + AXISDEVIATION_FILE ); + + if ( errors > 5 ) + error( "read_axisdev", "Exiting."); + + // try to re-generate the file + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_axisdev", "Generating file: %s\n", + AXISDEVIATION_FILE ); + + // try to open the file + axisdevfile_write.open( AXISDEVIATION_FILE ); + + // write short comment in the beginning + + axisdevfile_write << "# axis deviation file" << endl + << "# axis deviation for each mirror" << endl + << "# J C Gonzalez, 1998" << endl + << "#" << endl; + + // write data + + for ( i=0; i= VERBOSE_MINIMAL ) { + if ( errors > 0 ) + log( "read_axisdev", "Reading axis deviation file after %d fail%c\n", + errors, ((errors>1)?'s':' ') ); + else + log( "read_axisdev", "Reading axis deviation file . . .\n" ); + } + + //allocate memory for AxisDeviation table + + AxisDeviation = new float * [2]; + AxisDeviation[0] = new float[ct_NMirrors]; + AxisDeviation[1] = new float[ct_NMirrors]; + + // scan the whole file + + i = 0; + + while ( (! axisdevfile.eof()) && (i>>>>>>>>> MAGIC <<<<<<<<<< + if ( ct_Type == 1) { + + // allocate memory for the focals table + ct_Focal = new float[ct_NMirrors]; + + for (i=0; i>>>>>>>>> MAGIC <<<<<<<<<< + + // try to open the file + focalfile.open( FOCALS_FILE ); + + // if it is wrong or does not exist: + + while ( focalfile.bad() ) { + + ++errors; + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_focals", "Cannot open focals file: %s\n", + FOCALS_FILE ); + + if ( errors > 5 ) + error( "read_focals", "Exiting."); + + // try to re-generate the file + + if ( verbose >= VERBOSE_MINIMAL ) + log( "read_focals", "Generating file: %s\n", + FOCALS_FILE ); + + // try to open the file + focalfile_write.open( FOCALS_FILE ); + + // write short comment in the beginning + + focalfile_write << "# focals file" << endl + << "# focals for each mirror in the frame" << endl + << "# J C Gonzalez, 1998" << endl; + + // write data + + rnormal( NormalRandomNumbers, ct_NMirrors+1 ); + + for (i=0; i= VERBOSE_MINIMAL ) { + + if ( errors > 0 ) + log( "read_focals", "Reading focals file after %d fail%c\n", + errors, ((errors>1) ? 's' : ' ') ); + + else + + log( "read_focals", "Reading focals file . . .\n" ); + } + + // scan the whole file + + i = 0; + + while ( (! focalfile.eof()) && (iEtotal, (int)evth->EvtNumber); + + // CRITERIA + energy = evth->Etotal; + + if ( (energyEup) ) { + if (verbose) + log("REJECTED", "\n"); + reject = TRUE; + } else { + if (verbose) + log("ACCEPTED", "\n"); + reject = FALSE; + } + + if ( !reject ) { + + // update no. of cer. file + sprintf(cername, "%s/cer%06d", TMP_STDIN_DIR, ncerf+1); + + cerfile.open( cername ); + + if ( verbose ) + log("get_stdin_files", " Writing event %d, will be %s\n", + (int)evth->EvtNumber, cername); + + // write it + cerfile.write( (char*)buffer, BUFFER_LENGTH * sizeof(float) ); + + } + + // read blocks till new EVTH is found + // (now must be at the beginning of a block) + + if (verbose) + log("get_stdin_files", " Reading photons...\n"); + + bytes = 0; + while ( TRUE ) { + + // read photon + cin.read( (char*)buffer, 7 * sizeof(float) ); + + // the last photon is reached when the EVTH is read, but + // only when bytes % SIZE_OF_BLOCK == 0 + // SIZE_OF_BLOCK = 22932 + // with this, we only execute strstr once each 819 photons + + // are we in a new block? + + if ((bytes % SIZE_OF_BLOCK) == 0) { + // is this the last photon? + if ( strstr(chkbuffer, EVTH) == chkbuffer ) + break; + } + + // if not, write it to the cer-file + if ( !reject ) + cerfile.write( (char*)buffer, 7 * sizeof(float) ); + + bytes += 7 * sizeof(float); + + } + + if (verbose) + log("get_stdin_files", " Finished: %d bytes\n", bytes); + + if ( !reject ) { + // close cer-file + cerfile.close(); + } + + // at this point, we must be reading sta-file + + // read the rest of the sta-file + if (verbose) + log("get_stdin_files", " Reading sta-data...\n"); + cin.read( (char*)buffer + 7*sizeof(float), 5956 - 7*sizeof(float)); + + if ( !reject ) { + + // log + if (verbose) + log("get_stdin_files", " Saving sta-file...\n"); + + // open file + sprintf(staname, "%s/sta%06d", TMP_STDIN_DIR, ncerf+1); + stafile.open( staname ); + + // write it + stafile.write( (char*)buffer, 5956 ); + + // close sta-file + stafile.close(); + + // get out from here, return to main program + return( TRUE ); + + } + + } + + return( TRUE ); + +} +//!@} + + +//!--------------------------------------------------------------------- +// @name get_new_ct_pointing +// +// @desc returns new random position of CT around a given one +// +// @var theta Theta (ZA) of the shower +// @var phi Phi (AZ) of the shower +// @var range Maximum allowed distance between CT and shower axis +// @var newtheta Resulting Theta of the CT +// @var newphi Resulting Phi of the CT +// +// @return Angular distance between orig. and new directions +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function + +//!@{ +float +get_new_ct_pointing(float theta, float phi, float range, + float *newtheta, float *newphi) +{ + float distance; + float it, ip; + float sin_theta, cos_theta; + float sin_newtheta, cos_newtheta; + float sin_iphi, cos_iphi; + float iphi; + + // for the moment, we only simulate an uniform distribution, + // since our theta distribution in the generation of events is + // already uniform for hadrons, which are the main targets for + // using this option + + it = RandomNumber * range; + ip = RandomNumber * 2.0 * M_PI; + + sin_theta = sin(theta); + cos_theta = cos(theta); + + cos_newtheta = cos_theta*cos(it) + sin_theta*sin(it)*cos(ip); + *newtheta = acos( cos_newtheta ); + sin_newtheta = sin( *newtheta ); + + sin_iphi = sin(it)*sin(ip) / sin_newtheta; + cos_iphi = (( cos(it) - cos_newtheta * cos_theta ) / + ( sin_newtheta * sin_theta )); + + iphi = atan2( sin_iphi, cos_iphi ); + + *newphi = phi + iphi; + + return( it ); +} +//!@} + + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ +// +// $Log: not supported by cvs2svn $ +// Revision 1.18 1999/10/05 11:11:12 gonzalez +// Sep.1999 +// +// Revision 1.17 1999/03/24 16:33:01 gonzalez +// REFLECTOR 1.1: Release +// +// Revision 1.16 1999/01/21 16:03:44 gonzalez +// Only small modifications +// +// Revision 1.15 1999/01/19 18:07:16 gonzalez +// Bugs in STDIN-STDOUT version corrected. +// +// Revision 1.14 1999/01/14 17:35:47 gonzalez +// Both reading from STDIN (data_from_stdin) and +// writing to STDOUT (data_to_STDOUT) working. +// +// Revision 1.13 1999/01/13 12:41:12 gonzalez +// THIS IS A WORKING (last?) VERSION +// +// Revision 1.12 1998/12/17 16:26:09 gonzalez +// ************************************************* +// ************************************************* +// WARNING:: Version 1.11 is completly wrong!! +// ************************************************* +// ************************************************* +// +// Revision 1.10 1998/12/15 10:47:22 gonzalez +// RELEASE-1.0 +// +// Revision 1.9 1998/11/25 16:30:47 gonzalez +// Commit after inclusion of 'Blocking' +// +//!@} + +//=EOF Index: trunk/MagicSoft/Simulation/Detector/Reflector/reflector.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/Reflector/reflector.h (revision 292) +++ trunk/MagicSoft/Simulation/Detector/Reflector/reflector.h (revision 292) @@ -0,0 +1,277 @@ +//=////////////////////////////////////////////////////////////////////// +//= +//= reflector +//= +//= @file reflector.h +//= @desc Header file +//= @author J C Gonzalez +//= @email gonzalez@mppmu.mpg.de +//= @date Thu May 7 16:24:22 1998 +//= +//=---------------------------------------------------------------------- +//= +//= Created: Thu May 7 16:24:22 1998 +//= Author: Jose Carlos Gonzalez +//= Purpose: Program for reflector simulation +//= Notes: +//= +//=---------------------------------------------------------------------- +//= +//= $RCSfile: reflector.h,v $ +//= $Revision: 1.1.1.1 $ +//= $Author: harald $ +//= $Date: 1999-10-29 07:00:34 $ +//= +//=////////////////////////////////////////////////////////////////////// + +// @T \newpage + +//!@section Source code of |reflector.h|. + +/*!@" + + This section shows the include file |reflector.h| + +@"*/ + +//!@{ + +#ifndef _this_ +#define _this_ readparam +#endif + +//!@} + +//=----------------------------------------------------------- +//!@subsection Include files. + +//!@{ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "jcmacros.h" +#include "jcdebug.h" + +#include "reflector-v.h" + +#include "readparam.h" +#include "atm.h" + +#include "lagrange.h" + +#include "COREventHeader.hxx" +#include "CORParticle.hxx" +#include "CORStatfile.hxx" +#include "MCEventHeader.hxx" +#include "MCCphoton.hxx" + +// command line options available +#define COMMAND_LINE_OPTIONS "f:h" + +// This is C++, but RANLIB routines are written in pure ANSI C. +// In order to read easily these routines, we must include +// the following directive +extern "C" { + +#include "ranlib.h" + +} + +//!@} + +// @subsection Macro-definitions, and constants + +//!@{ + +// random numbers +#define RandomNumber ranf() + +// Speed of Light in vacuum, in m/s +const float Speed_of_Light_vacuum = 299792458.0; // EXACT!! +const float Speed_of_Light_air = Speed_of_Light_vacuum / 1.000293; + +// Speed of Light in vacuum, in cm/ns +const float Speed_of_Light_vacuum_cmns = Speed_of_Light_vacuum / 1.0e7; +const float Speed_of_Light_air_cmns = Speed_of_Light_air / 1.0e7; + +// now we define the list CT_ITEM_LIST of possible items in +// the CT definition file. + +#define CT_ITEM_LIST /* LIST OF ITEMS IN THE CT DEFINITION FILE */ \ +T(type), /* type of definition file */ \ +T(focal_distance), /* focal distance */ \ +T(focal_std), /* std(focal distance) */ \ +T(point_spread), /* point spread */ \ +T(point_std), /* std(point spread) */ \ +T(adjustment_dev), /* std of adjustment deviation */ \ +T(black_spot), /* radius of the black spot in center of mirrors */ \ +T(n_mirrors), /* number of mirrors */ \ +T(r_mirror), /* radius of one mirror */ \ +T(camera_width), /* camera width */ \ +T(n_pixels), /* number of pixels in the camera */ \ +T(pixel_width), /* pixel width */ \ +T(define_mirrors) /* this entry is followed by the def. of pixels */ + +#define T(x) x // define T() as the name as it is + + enum CT_ITEM_TYPE { + CT_ITEM_LIST + }; + +#undef T + +#define T(x) #x // define T() as the string of x + +const char *const CT_ITEM_NAMES[] = { + CT_ITEM_LIST +}; + +#undef T + +// TYPE=0 (CT1) +// i s rho theta x y z thetan phin xn yn zn +// +// i : number of the mirror +// s : arc length [cm] +// rho : polar rho of the position of the center of the mirror [cm] +// theta : polar angle of the position of the center of the mirror [cm] +// x : x coordinate of the center of the mirror [cm] +// y : y coordinate of the center of the mirror [cm] +// z : z coordinate of the center of the mirror [cm] +// thetan : polar theta angle of the direction where the mirror points to +// phin : polar phi angle of the direction where the mirror points to +// xn : xn coordinate of the normal vector in the center (normalized) +// yn : yn coordinate of the normal vector in the center (normalized) +// zn : zn coordinate of the normal vector in the center (normalized) +// +// TYPE=1 (MAGIC) +// i f sx sy x y z thetan phin +// +// i : number of the mirror +// f : focal distance of that mirror +// sx : curvilinear coordinate of mirror's center in X[cm] +// sy : curvilinear coordinate of mirror's center in X[cm] +// x : x coordinate of the center of the mirror [cm] +// y : y coordinate of the center of the mirror [cm] +// z : z coordinate of the center of the mirror [cm] +// thetan : polar theta angle of the direction where the mirror points to +// phin : polar phi angle of the direction where the mirror points to +// xn : xn coordinate of the normal vector in the center (normalized) +// yn : yn coordinate of the normal vector in the center (normalized) +// zn : zn coordinate of the normal vector in the center (normalized) + +#define CT_I 0 + +#define CT_S 1 +#define CT_RHO 2 +#define CT_THETA 3 + +#define CT_FOCAL 1 +#define CT_SX 2 +#define CT_SY 3 + +#define CT_X 4 +#define CT_Y 5 +#define CT_Z 6 +#define CT_THETAN 7 +#define CT_PHIN 8 +#define CT_XC 9 +#define CT_YC 10 +#define CT_ZC 11 + +#define CT_NDATA 12 + +// pre-defined filenames / values + +#define REFLECTIVITY_FILE "reflectivity.dat" +#define FOCALS_FILE "focals.dat" +#define AXISDEVIATION_FILE "axisdev.dat" + +// values for read-ot from STDIN + +#define BUFFER_LENGTH 5733 +#define SIZE_OF_BLOCK 22932 +static char EVTH[] = "EVTH"; + +// maximum deviation from the original direction with +// "random_pointing" option +// +// now is 6 degrees = 0.104719755119660 radians + +#define RANDOM_POINTING_MAX_SEPARATION 0.104719755119660 + +//!@} + +//=----------------------------------------------------------- +//!@subsection Prototypes of functions. + +//!@{ + +//++ +// prototypes +//-- + +void present(void); +void usage(void); +void clean(void); +void log(const char *funct, char *fmt, ...); +void logerr(const char *funct, char *fmt, ...); +void error(const char *funct, char *fmt, ...); +void makeOmega(float theta, float phi); +void makeOmegaI(float theta, float phi); +void applyMxV(float M[3][3], float *V, float *Vp); +void read_ct_file(void); +void read_reflectivity(void); +void read_axisdev(void); +void read_focals(void); +void rnormal(double *r, int n); +float dist_r_P(float a, float b, float c, + float u, float v, float w, + float x, float y, float z); +void write_mark(ofstream &os, char mark=0xff, int length=1); +int isA( char * s1, const char * flag ); +float Curv2Lin(float s); +float Lin2Curv(float x); +int get_stdin_files(int ncerf, float El=0., float Eu=1000000., + int flag=FALSE); +float get_new_ct_pointing(float theta, float phi, float range, + float *newtheta, float *newphi); + +//!@} + +//=------------------------------------------------------------ +//!@subsection Log of this file. + +//!@{ + +/* + * $Log: not supported by cvs2svn $ + * Revision 1.11 1999/10/05 11:11:12 gonzalez + * Sep.1999 + * + * Revision 1.10 1999/03/24 16:33:03 gonzalez + * REFLECTOR 1.1: Release + * + * Revision 1.9 1999/01/21 16:03:53 gonzalez + * Only small modifications + * + * Revision 1.8 1999/01/19 18:07:17 gonzalez + * Bugs in STDIN-STDOUT version corrected. + * + * Revision 1.7 1999/01/14 17:35:44 gonzalez + * Both reading from STDIN (data_from_stdin) and + * writing to STDOUT (data_to_STDOUT) working. + * + */ + +//!@} +//=EOF