Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile (revision 725) @@ -0,0 +1,211 @@ +################################################################## +# +# 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: 2001-04-09 08:59:01 $ +# +################################################################## +# @maintitle + +# @code + +INCLUDEMK = config.mk.${OSTYPE} +include ${INCLUDEMK} + +# @endcode + +# @code + +# common flags +INCLUDES = -I/usr/include -I/usr/include/cxx + +#CERNLIBDIR = ${CERNDIR}/pro/lib/ +#CERNLIB = -L${CERNLIBDIR} -lgraflib -lgrafX11 -lpacklib -lkernlib -lpawlib + +RANLIB = -L${RANLIBDIR} -lranlib + +# special flags + +osf_FORLIBS = -L/usr/lib -lUfor -lFutil -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} +CFLAGS = -std1 -warnprotos -msg_enable level6 -msg_disable returnchecks +FFLAGS = ${CXXFLAGS} +LIBS = ${RANLIB} ${FORLIBS} +LIBS = ${CERNLIB} ${RANLIB} ${FORLIBS} + +#------------------------------------------------------------------------------ + +#.SILENT: + +.SUFFIXES: .c .cxx .C .c++ .h .hxx .H .h++ .o .so .f + +SRCS = \ + attenu.f \ + diag.c \ + init.c \ + parms.c \ + geometry.c \ + atm.c \ + ph2cph.c \ + header.c \ + reflector.c + +HEADERS = \ + atm.h \ + diag.h \ + geometry.h \ + header.h \ + init.h \ + lagrange.h \ + parms.h \ + version.h + +OBJS = \ + attenu.o \ + diag.o \ + init.o \ + parms.o \ + geometry.o \ + atm.o \ + ph2cph.o \ + header.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 Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c (revision 725) @@ -0,0 +1,58 @@ +#include +#include +#include +#include "diag.h" +#include "atm.h" +#include "init.h" + +/* random numbers */ +#define RandomNumber ranf() + +/* Local declarations */ +static int atmModel=0; /* current atm. model */ + +/* Function declarations */ +static float atm(float wavelength, float height, float theta); +void SetAtmModel(char *model); +int absorption(float wlen, float height, float theta); +extern void attenu_(float *, float *, float *, float *); /* in Fortran */ +extern float ranf(void); + +void SetAtmModel(char *model) +{ + while (strcmp(model, AtmModelNames[atmModel])) + if (++atmModel == sizeof(AtmModelNames)/sizeof(AtmModelNames[0])) + { atmModel = 0; + Error(ATM__NFND__ERR, model); + break; } + + Log(ATM__SET___LOG, AtmModelNames[atmModel]); +} /* end of SetAtmModel */ + +static float atm(float wavelength, float height, float theta) +{ float transmittance = 1.0; /* final atm transmittance (ret. value) */ + + switch(atmModel) + { case ATM_NOATMOSPHERE: /* no atm at all: transmittance = 100% */ + break; + case ATM_90PERCENT: /* atm. with transmittance = 90% */ + transmittance = 0.9; + break; + case ATM_ISOTHERMAL: /* isothermal approximation */ + /********************/ + break; + case ATM_CORSIKA: /* atmosphere as defined in CORSIKA */ + attenu_(&wavelength, &height, &theta, &transmittance); + break; + } /* end of atm switch */ + + return transmittance; +} /* end of atm */ + +int absorption(float wlen, float height, float theta) +{ int ret = 0; /* 0: passed, 1: absorbed */ + + if ( RandomNumber > atm(wlen, height, theta)) return 1; + + return ret; +} /* end of absorption */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h (revision 725) @@ -0,0 +1,25 @@ +#ifndef __RFL_ATM__ +#define __RFL_ATM__ + +#define ITEM_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 { ITEM_LIST }; +#undef T + +#define T(x) #x /* define T() as the string of x */ + const char *AtmModelNames[] = { ITEM_LIST }; +#undef T +#undef ITEM_LIST + +/* Strings */ +#define ATM__NFND__ERR /* model name */ \ + " *** Atm model \"%s\" not found.\n" +#define ATM__SET___LOG /* model name */ \ + "Setting atm model \"%s\".\n" + +#endif Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f (revision 725) @@ -0,0 +1,471 @@ +********************************************************************* +* * +* 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 + + DOUBLE PRECISION Rcos2, Rsin2 + DOUBLE PRECISION x1, x2, x3, x4 + DOUBLE PRECISION e1, e2, e3, e4 + + REAL wavelength, height, theta, tr_atmos + real trr, trm, tro + REAL LONG(12) +c fs: define obervation level + double precision obslev + 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 / + +c-- Take same Earth radius as in CORSIKA + parameter (rt = 6371315.D2) + +c-- Scale-height for Exponential density profile + parameter (hscale = 7.4d5) + +*********************************************************************** +* +* SCATTERING PARAMETERS FOR RAYLEIGH: +* +* MEAN FREE PATH FOR SCATTERING RAYLEIGH (g/cm^2) + PARAMETER (XR=2.970D3) +*********************************************************************** + +c-- Observation level at La Palma + parameter (obslev = 2200.d2) + + T_Ray = 1.0 + T_Mie = 1.0 + T_Oz = 1.0 + +c-- Height calculated using an obslev = H(LaPalma) <> 0 +c H = -RT + SQRT(RT**2 + (height/COS(theta))**2 + +c + (2.0D0*RT*height)) + h = -rt + sqrt((rt+obslev)**2 + + + ((height-obslev)/cos(theta))**2 + + + (2.0d0*(rt+obslev)*(height-obslev))) + + ROW = AINT(((H+1.)/1.0E5)) + +*********************************************************************** +* +* LARGE ZENITH ANGLE FACTOR (AIR MASS): + +c fs : air mass factor = path_lenght(za) / path_lenght(vertical) +c at point of emission of cherenkov photon +c => pure geometric correction on +c absorption lenght measured at vertical height [km^-1] +c-- +c ma = (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)- +c + (RT*COS(theta))) +c-- +c mb = (1.d0/(h-obslev)) * +c & ( -(rt+obslev)*cos(theta) +c & +sqrt( (rt+h)**2 - ((rt+obslev)*sin(theta))**2) ) +c-- + +c-- Air mass "m" calcualted using a one-exponential density +c-- profile for the atmosphere, rho = rho_0 exp(-h/Hs) +c-- with Hs = hscale = 7.4 km + + Rcos2 = rt * cos(theta)**2 + Rsin2 = rt * sin(theta)**2 + + x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale)) + x2 = sqrt((2. * h + Rcos2) / (2. * hscale)) + x3 = sqrt((2. * obslev + rt) / (2. * hscale)) + x4 = sqrt((2. * h + rt) / (2. * hscale)) + + e1 = derfc(x1) + e2 = derfc(x2) + e3 = derfc(x3) + e4 = derfc(x4) + + m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4)) + +********************************************************************** +* +* 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/ReflectorII/config.mk.linux-gnu =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/config.mk.linux-gnu (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/config.mk.linux-gnu (revision 725) @@ -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.linux-gnu,v $ +# $Revision: 1.1.1.1 $ +# $Author: harald $ +# $Date: 2001-04-09 08:59:10 $ +################################################################## +# @maintitle + +# @code + +# program + +PROGRAM = camera + +# compilers & tools + +CC = gcc +CXX = g++ +F77 = g77 + +DOCUM = ../sus/sus +RATE = ../../comstat/comstat -g -p + +# includes + +INCLUDE = ../include-GENERAL +INCLUDE_COR = ../include-CORSIKA +INCLUDE_MC = ../include-MC +INCLUDE_REFL = . + +OPTIM = +DEBUG = -g + +# libraries + +RANLIBDIR = ../lib +CERNDIR = /CERN + +# system + +SYSTEM = linux + +# @endcode +##EOF Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/config.mk.osf1 =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/config.mk.osf1 (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/config.mk.osf1 (revision 725) @@ -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.osf1,v $ +# $Revision: 1.1.1.1 $ +# $Author: harald $ +# $Date: 2001-04-09 08:59:10 $ +################################################################## +# @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/ReflectorII/diag.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/diag.c (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/diag.c (revision 725) @@ -0,0 +1,94 @@ +#include +#include +#include +#include "version.h" +#include "diag.h" + +/* To display errors on stdout, comment line below */ +#define ERR_TO_STDERR + +#ifdef ERR_TO_STDERR + #define STDERR stderr +#elif + #define STDERR stdout +#endif + +enum VerboseLevel /* verbose levels */ +{ VERBOSE_QUIET, + VERBOSE_MINIMAL, + VERBOSE_NORMAL, + VERBOSE_MAXIMAL }; + +#define VERBOSE_DEFAULT VERBOSE_QUIET + +static int verbose = VERBOSE_DEFAULT; + +void SetVerbose(int vlevel) +{ verbose = (vlevel < VERBOSE_QUIET) ? VERBOSE_QUIET + : (vlevel > VERBOSE_MAXIMAL) ? VERBOSE_MAXIMAL + : vlevel; +} /* end of SetVerbose */ + +void FatalError(const char *fatal_string, ...) +{ va_list args; + extern void clean(void); + + /* Display signature */ + fprintf(STDERR, "\a\n *** %s %s: ", QUOTE(PROGRAM), QUOTE(VERSION)); + + /* Display error message */ + va_start(args, fatal_string); + vfprintf(STDERR, fatal_string, args); + va_end(args); + + /* Byebye message */ + fputs(" *** Exiting Reflector.\n\n", STDERR); + + /* Abort and exit */ + clean(); + exit(1); +} /* end of FatalError */ + +void Error(const char *error_string, ...) +{ va_list args; + + /* Display signature */ + fprintf(STDERR, "\a\n *** %s %s: ", QUOTE(PROGRAM), QUOTE(VERSION)); + + /* Display error message */ + va_start(args, error_string); + vfprintf(STDERR, error_string, args); + va_end(args); + + /* Skip line */ + fputc('\n', STDERR); +} /* end of Error */ + +void Message(const char *message_string, ...) +{ va_list args; + + /* Display message */ + va_start(args, message_string); + vprintf(message_string, args); + va_end(args); +} /* end of Message */ + +void Log(const char *log_string, ...) +{ va_list args; + + /* Display message */ + va_start(args, log_string); + if (verbose >= VERBOSE_MINIMAL) + vprintf(log_string, args); + va_end(args); +} /* end of Log */ + +void Debug(const char *log_string, ...) +{ va_list args; + + /* Display message */ + va_start(args, log_string); + if (verbose == VERBOSE_MAXIMAL) + vprintf(log_string, args); + va_end(args); +} /* end of Debug */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/diag.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/diag.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/diag.h (revision 725) @@ -0,0 +1,10 @@ +#ifndef __RFL_DIAG__ +#define __RFL_DIAG__ + +void FatalError(const char *fatal_string, ...); +void Error(const char *error_string, ...); +void Message(const char *message_string, ...); +void Log(const char *log_string, ...); +void Debug(const char *log_string, ...); + +#endif Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.c (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.c (revision 725) @@ -0,0 +1,213 @@ +#include +#include +#include +#include "diag.h" +#include "geometry.h" +#include "init.h" + +extern char line[]; /* parsing buf. (init) */ + +float ct_Focal_mean; /* focal dist. (mean) (cm) */ +float ct_Focal_std; /* focal dist. (std) (cm) */ +float ct_PSpread_mean; /* pt. spread fn. (mean) (cm) */ +float ct_PSpread_std; /* pt. spread fn. (std) (cm) */ +float ct_Adjustment_std; /* adjustment dev. (std) (cm) */ +float ct_BlackSpot_rad; /* black spot radius (cm) */ +float ct_RMirror; /* rad. of single mirror (cm) */ +int ct_NMirrors=0; /* number of mirrors */ +float ct_CameraWidth; /* camera width (cm) */ +int ct_NPixels; /* number of pixels */ +float ct_PixelWidth; /* pixel width (cm) */ + +mirror *ct_data=NULL; /* ptr to mirror data */ +FILE *ct_BinaryData=NULL; /* binary data for mirrors */ +char ct_BinaryName[128]; /* binary data filename */ + +int nReflectivity=0; /* elements in refl. table */ +float *Reflectivity[2]; /* reflectivity table */ +float *AxisDeviation[2]; /* axis deviation table */ +float *ct_Focal; /* focal length table */ + +/* Prototypes */ +static void ReadMirrorTable(FILE *geofile); +static void ReadReflectivity(char *datname); +static void ReadAxisDev(char *datname); +static void ReadFocals(void); + +static void ReadMirrorTable(FILE *geofile) +{ int i; /* Mirror index */ + + if ((ct_data=(mirror *)malloc(sizeof(mirror)*ct_NMirrors)) == NULL) + FatalError(MIRR_ALLOC_FTL, ct_NMirrors); + Log(MIRR_ALLOC_LOG, ct_NMirrors); + Log(MIRR_TABLE_LOG); + if (ct_BinaryData) + { Log(BINF_OPEN__LOG, ct_BinaryName); + fread(ct_data, sizeof(mirror), ct_NMirrors, ct_BinaryData); + fclose(ct_BinaryData); } + else + { /* read ASCII data */ + Log(READ_ASCII_LOG); + for (i=0; i= nReflectivity) break; }} + fclose(datfile); + + nReflectivity = current; +} /* end of ReadReflectivity */ + +static void ReadAxisDev(char *datname) +{ FILE *datfile = fopen(datname, "r"); + int current = 0; + + if (datfile == NULL) + FatalError(AXIS_ERROR_FTL, datname); + if ((AxisDeviation[0]= + (float *) malloc(sizeof(float) * ct_NMirrors)) == NULL + || (AxisDeviation[1]= + (float *) malloc(sizeof(float) * ct_NMirrors)) == NULL) + FatalError(AXIS_ALLOC_FTL, ct_NMirrors); + while (fgets(line, LINE_MAX_LENGTH, datfile)) + { if (line[0] == '#') continue; + if (2==sscanf(line, "%f %f", + &AxisDeviation[0][current], &AxisDeviation[1][current])); + { current++; + if (current >= ct_NMirrors) break; }} + fclose(datfile); + + if (current != ct_NMirrors) + FatalError(AXIS_FEW___FTL, current, ct_NMirrors); +} /* end of ReadAxisDev */ + +static void ReadFocals(void) +{ int loop; + + if ((ct_Focal = (float *) malloc(sizeof(float) * ct_NMirrors)) == NULL) + FatalError(FOCL_FEW___FTL, ct_NMirrors); + + for (loop=0; loopparm val. */ + + /* Initialise arrays */ + Reflectivity[0] = AxisDeviation[0] = NULL; + + do + { switch(ParseLine(geofile, ctparms, ARRAY_SZ(ctparms), &value_ptr)) + { case type: + if (1 != atoi(value_ptr)) + FatalError(TYPE_ERROR_FTL); + break; + case focal_distance: + Log(LOG__FLOAT_LOG, "focal distance (average)", + ct_Focal_mean = (float) atof(value_ptr)); + break; + case focal_std: + Log(LOG__FLOAT_LOG, "focal distance (std. dev.)", + ct_Focal_std = (float) atof(value_ptr)); + break; + case point_spread: + Log(LOG__FLOAT_LOG, "point spread fn. (average)", + ct_PSpread_mean = (float) atof(value_ptr)); + break; + case point_std: + Log(LOG__FLOAT_LOG, "point spread fn. (std. dev.)", + ct_PSpread_std = (float) atof(value_ptr)); + break; + case adjustment_dev: + Log(LOG__FLOAT_LOG, "adjustment dev. (std. dev.)", + ct_Adjustment_std = (float) atof(value_ptr)); + break; + case black_spot: + Log(LOG__FLOAT_LOG, "radius of black spot (cm)", + ct_BlackSpot_rad = (float) atof(value_ptr)); + break; + case r_mirror: + Log(LOG__FLOAT_LOG, "single mirror radius (cm)", + ct_RMirror = (float) atof(value_ptr)); + break; + case n_mirrors: + Log(LOG__INT___LOG, "number of mirrors", + ct_NMirrors = atoi(value_ptr)); + break; + case camera_width: + Log(LOG__FLOAT_LOG, "camera width (cm)", + ct_CameraWidth = (float) atof(value_ptr)); + break; + case n_pixels: + Log(LOG__INT___LOG, "number of pixels", + ct_NPixels = atoi(value_ptr)); + break; + case pixel_width: + Log(LOG__FLOAT_LOG, "pixel width (cm)", + ct_PixelWidth = (float) atof(value_ptr)); + break; + case refl_file: + ReadReflectivity(value_ptr); + break; + case axisdev_file: + ReadAxisDev(value_ptr); + break; + case define_mirrors: + if (ct_NMirrors) ReadMirrorTable(geofile); + else FatalError(MIRR_NSPEC_FTL); + switch_end = TRUE; + break; + default: switch_end = TRUE; + break; }} + while (!switch_end); + fclose(geofile); + + if (Reflectivity[0] == NULL) ReadReflectivity(REFLECTIVITY_FILE); + if (AxisDeviation[0]== NULL) ReadAxisDev(AXISDEVIATION_FILE); + ReadFocals(); +} /* end of GeometrySwitch */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.h (revision 725) @@ -0,0 +1,86 @@ +#ifndef __RFL_GEOMETRY__ +#define __RFL_GEOMETRY__ + +#define 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(refl_file), /* path of file containing refl. data */ \ +T(axisdev_file), /* path of file containing axis dev. data */ \ +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 { ITEM_LIST }; +#undef T + +#define T(x) #x /* define T() as the string of x */ + const char *ctparms[] = { ITEM_LIST }; +#undef T +#undef ITEM_LIST + +/* Strings */ +#define TYPE_ERROR_FTL /* no parms */ \ + "This version of Reflector handles only MAGIC.\n" \ + " *** Please change \"ct_type\" in geometry file accordingly.\n" +#define MIRR_NSPEC_FTL /* no parms */ \ + "The number of mirrors is needed for further processing.\n" \ + " *** Please specify it after the \"n_mirrors\" directive.\n" +#define LOG__FLOAT_LOG /* floatname, value */ \ + "<%s> set at %.2f.\n" +#define LOG__INT___LOG /* intname, value */ \ + "<%s> set at %d.\n" +#define MIRR_TABLE_LOG /* no parms */ \ + "Start reading mirror data.\n" +#define MIRR_ALLOC_FTL /* nr of mirrors */ \ + "Cannot allocate enough memory space for %d mirrors.\n" +#define MIRR_ALLOC_LOG /* nr of mirrors */ \ + "Allocated memory for %d mirrors.\n" +#define BINF_OPEN__LOG /* bin. data filename */ \ + "Located and opened file \"%s\" containing mirror data.\n" +#define READ_ASCII_LOG /* no parms */ \ + "Start reading ASCII data for mirror nr:\n" +#define MIRR_FEW___FTL /* mirrors read */ \ + "Not enough mirror data: only %d mirrors read.\n" +#define BINF_ERROR_LOG /* bin. data filename */ \ + "Cannot write mirror data on binary file:\n %s\n" \ + " *** Skipping operation and proceeding.\n" +#define BINF_WRITE_LOG /* bin. data filename */ \ + "Binary data written in file \"%s\".\n" +#define RFLF_ERROR_FTL /* reflectivity fname */ \ + "Cannot find file \"%s\" containing data on reflectivity.\n" +#define REFL_ALLOC_FTL /* nr of refl items */ \ + "Cannot allocate enough memory space for %d refl. data.\n" +#define AXIS_ERROR_FTL /* axis dev fname */ \ + "Cannot find file \"%s\" containing data on axis dev.\n" +#define AXIS_ALLOC_FTL /* nr of a/d items */ \ + "Cannot allocate enough memory space for %d axis dev. data.\n" +#define AXIS_FEW___FTL /* read, total items */ \ + "Found only %d out of %d axis dev. data.\n" +#define FOCL_FEW___FTL /* nr of focal items */ \ + "Cannot allocate enough memory space for %d focal data.\n" + +/* pre-defined filenames / values */ + +#define REFLECTIVITY_FILE "../Data/reflectivity.dat" +#define AXISDEVIATION_FILE "../Data/axisdev.dat" + +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 + +#endif Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/header.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/header.c (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/header.c (revision 725) @@ -0,0 +1,43 @@ +#include +#include "header.h" + +/* Need to be defined the procedures to set the values in + float TimeFirst; + float TimeLast; */ + +static RflHeader rhead; RflHeader *rheadp = &rhead; +static CerHeader chead; CerHeader *cheadp = &chead; + +void TranslateHeader(RflHeader *r, CerHeader *c) +{ + r->EvtNumber = c->EvtNumber; + r->PrimaryID = c->PrimaryID; + r->Etotal = c->Etotal; + r->Thick0 = c->Thick0; + r->FirstTarget = c->FirstTarget; + r->zFirstInt = c->zFirstInt; + r->Theta = c->Theta; + r->Phi = c->Phi; + r->NumRndSeq = c->NumRndSeq; + r->RunNumber = c->RunNumber; + r->DateRun = c->DateRun; + r->VersionPGM = c->VersionPGM; + r->NumObsLev = c->NumObsLev; + r->SlopeSpec = c->SlopeSpec; + r->ELowLim = c->ELowLim; + r->EUppLim = c->EUppLim; + r->ThetaMin = c->ThetaMin; + r->ThetaMax = c->ThetaMax; + r->PhiMin = c->PhiMin; + r->PhiMax = c->PhiMax; + r->CWaveLower = c->CWaveLower; + r->CWaveUpper = c->CWaveUpper; + + memcpy(r->p, c->p, 3*sizeof(float)); + memcpy(r->RndData, c->RndData, 30*sizeof(float)); + memcpy(r->HeightLev, c->HeightLev, 10*sizeof(float)); + memcpy(r->CorePos, c->CorePos, 40*sizeof(float)); + + r->deviationPhi = r->deviationTheta = r->Trigger = 0.f; + +} /* end of TranslateHeader */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h (revision 725) @@ -0,0 +1,132 @@ +#ifndef __RFL_HEADER__ +#define __RFL_HEADER__ + +typedef struct /* EVTH from cerfile */ +{ char EVTH[4]; + float EvtNumber; + float PrimaryID; + float Etotal; + float Thick0; + float FirstTarget; + float zFirstInt; + float p[3]; + float Theta; + float Phi; + + float NumRndSeq; + float RndData[10][3]; + + float RunNumber; + float DateRun; + float VersionPGM; + + float NumObsLev; + float HeightLev[10]; + + float SlopeSpec; + float ELowLim; + float EUppLim; + + float Ecutoffh; + float Ecutoffm; + float Ecutoffe; + float Ecutoffg; + + float NFLAIN; + float NFLDIF; + float NFLPI0; + float NFLPIF; + float NFLCHE; + float NFRAGM; + + float Bx; + float By; + + float EGS4yn; + float NKGyn; + float GHEISHAyn; + float VENUSyn; + float CERENKOVyn; + float NEUTRINOyn; + float HORIZONTyn; + float COMPUTER; + + float ThetaMin; + float ThetaMax; + float PhiMin; + float PhiMax; + + float CBunchSize; + float CDetInX,CDetInY; + float CSpacInX,CSpacInY; + float CLenInX,CLenInY; + float COutput; + + float AngleNorthX; + float MuonInfo; + + float StepLength; + float CWaveLower; + float CWaveUpper; + float Multipl; + float CorePos[2][20]; + + float dmmy1; + float SpinTheta; + float SpinPhi; + float dmmy2[132]; +} CerHeader; + +typedef struct /* EVTH from rflfile */ +{ float EvtNumber; + float PrimaryID; + float Etotal; + float Thick0; + float FirstTarget; + float zFirstInt; + float p[3]; + float Theta; + float Phi; + + float NumRndSeq; + float RndData[10][3]; + + float RunNumber; + float DateRun; + float VersionPGM; + + float NumObsLev; + float HeightLev[10]; + + float SlopeSpec; + float ELowLim; + float EUppLim; + + float ThetaMin; + float ThetaMax; + float PhiMin; + float PhiMax; + + float CWaveLower; + float CWaveUpper; + float CorePos[2][20]; + float TimeFirst; + float TimeLast; + + float deviationPhi; + float deviationTheta; + + float Trigger; + + float CORSIKAPhs; /* Original photons written by Corsika */ + float AtmAbsPhs; /* Photons absorbed by the atmosphere */ + float MirrAbsPhs; /* Photons absorbed by the mirror */ + float OutOfMirrPhs; /* Photons outside the mirror */ + float BlackSpotPhs; /* Photons lost in the "black spot" */ + float OutOfChamPhs; /* Photons outside the chamber */ + float CPhotons; /* Photons reaching the chamber */ +} RflHeader; + +void TranslateHeader(RflHeader *r, CerHeader *c); + +#endif Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c (revision 725) @@ -0,0 +1,115 @@ +#include +#include +#include +#include +#include "version.h" +#include "diag.h" +#include "init.h" + +/* extern declarations */ +FILE *filelist = NULL, /* filelist ptr */ + *rflname = NULL; /* reflector (out)file */ +char cername[LINE_MAX_LENGTH]; /* current cername */ +long first_Event = 0, /* first proc. event */ + last_Event = 1000000, /* last proc. event */ + max_Events = 50000; /* max proc. events */ +float low_Ecut = 0.0f, /* lower Ecut (GeV) */ + high_Ecut = 100000.0f; /* upper Ecut (GeV) */ + +int is_Fixed_Target = FALSE; /* fixed target? */ +float fixed_Theta, /* zenith angle (rad) */ + fixed_Phi; /* azi (0N->90E) (rad) */ + +int is_Random_Pointing=FALSE; /* random pointing? */ +float Random_Pointing_MaxDist; /* in metres */ +int nRepeat_Random = 1; /* nr of times a sh. is reused */ + +long Seeds[2] = {3141592L, 2718182L}; + +/* Other declarations */ +char whites[] = " \v\t\n\r"; /* White chars def */ +char line[LINE_MAX_LENGTH]; /* parsing buffer */ + +/* Prototypes */ +static void CheckSignature(char *line); +int ParseLine(FILE *parfile, /* FILE with parms */ + const char *token_list[], /* array w/tokens */ + int tokens, /* nr of tokens */ + char **value_ptr); /* ptr to parm val. */ + +/* Checks signature. + * It separates tokens and compare the first + * token with PROGRAM and the second with VERSION. + * Both values are stored in "version.h". + * Return value is set to true if check fails. */ + +static void CheckSignature(char *line) +{ char *cp; /* To parse signature */ + + if ((cp=strtok(line, whites)) == NULL || strcmp(cp, QUOTE(PROGRAM)) || + (cp=strtok(NULL, whites)) == NULL || strcmp(cp, QUOTE(VERSION))) + FatalError(SIGN_ERROR_FTL, QUOTE(PROGRAM), QUOTE(VERSION)); +} /* end of CheckSignature */ + +int ParseLine(FILE *parfile, /* FILE with parms */ + const char *token_list[], /* array w/tokens */ + int tokens, /* nr of tokens */ + char **value_ptr) /* ptr to parm val. */ +{ int item = tokens; + char *cp; + + do + { if (fgets(line, LINE_MAX_LENGTH, parfile) == NULL) + break; + else if (line[0] == '#'); /* skip comments (start with '#') */ + else if (line[0] == '>') /* show user comments (start with '>') */ + { if ((cp=strchr(line, '\n')) != NULL) *cp = 0; + puts(line); } + else if ((cp = strtok(line, whites)) != NULL) + { /* check which param it is */ + for (item=0; item= tokens); + + return item; +} /* end of ParseLine */ + +FILE *geofile = NULL; /* geometry file */ + +void init(char *filename) +{ FILE *parfile = NULL; /* parms file */ + extern void ParmsSwitch(FILE *parfile); /* std parsing fn */ + extern void GeometrySwitch(FILE *geofile); /* geometry parsing fn */ + + /* Open file or read from stdin */ + if (filename) parfile = fopen("filename", "r"); + if (parfile == NULL) parfile = stdin; + fgets(line, LINE_MAX_LENGTH, parfile); + + /* Check signature */ + CheckSignature(line); + + ParmsSwitch(parfile); + if (geofile) GeometrySwitch(geofile); + else FatalError(GEOM_NSPEC_FTL); + + /* Alloc memory for CPhotons */ + if ((CPhotons = + (cphoton *) malloc(NR_OF_CPHOTONS * sizeof(cphoton))) == NULL) + FatalError(CPHS_ALLOC_FTL, NR_OF_CPHOTONS); + else Log(CPHS_ALLOC_LOG, NR_OF_CPHOTONS); + + /* Write signature and START_OF_RUN at the beginning of rflfile */ + /* '\0' for backward compatibility */ + fseek(rflfile, 0L, SEEK_SET); + fprintf(rflfile, "%s %s", QUOTE(PROGRAM), QUOTE(VERSION)); + fputc(0, rflfile); + fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); + + /* "Splashscreen" */ + Message(RFL__START_MSG, QUOTE(PROGRAM), QUOTE(VERSION)); +} /* end of init */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h (revision 725) @@ -0,0 +1,136 @@ +#ifndef __RFL_INIT__ +#define __RFL_INIT__ + +/* Constant definitions */ +#define LINE_MAX_LENGTH 128 /* parsing buffer sz. */ +#define NR_OF_CPHOTONS 65535 /* CPhotons sz. */ +#define PH_IN_DATABLOCK 39 /* Photons in datablk */ +#define SHOW_ME 500 /* How many evts among two logs */ + +/* Macros */ +#define ARRAY_SZ(A) (sizeof(A)/sizeof(A[0])) + +/* Tags inside reflector file */ +#define SIZE_OF_FLAGS 13 + +#define FLAG_START_OF_RUN "\nSTART---RUN\n" +#define FLAG_START_OF_EVENT "\nSTART-EVENT\n" +#define FLAG_END_OF_EVENT "\nEND---EVENT\n" +#define FLAG_END_OF_RUN "\nEND-----RUN\n" +#define FLAG_END_OF_FILE "\nEND----FILE\n" +#define FLAG_END_OF_STDIN "\nEND---STDIN\n" + +/* Prototypes */ +void init(char *filename); + +/* External declarations */ +extern FILE *filelist; /* filelist ptr */ +extern FILE *rflfile; /* reflector (out)file */ +extern FILE *cerfile; /* current cerfile */ +extern char cername[]; /* current cername */ +extern long first_Event; /* first proc. event */ +extern long last_Event; /* last proc. event */ +extern long max_Events; /* max proc. events */ +extern float low_Ecut; /* lower Ecut (GeV) */ +extern float high_Ecut; /* upper Ecut (GeV) */ + +extern int is_Fixed_Target; /* fixed target? */ +extern float fixed_Theta; /* zenith angle (rad) */ +extern float fixed_Phi; /* azi (0N->90E) (rad) */ + +extern int is_Random_Pointing; /* random pointing? */ +extern float Random_Pointing_MaxDist; /* in metres */ +extern int nRepeat_Random; /* nr of times a sh. is reused */ +extern long Seeds[2]; /* random seeds */ + +extern FILE *ct_BinaryData; /* binary data for mirrors */ +extern char ct_BinaryName[]; /* binary data filename */ +extern float ct_Focal_mean; /* focal dist. (mean) (cm) */ +extern float ct_Focal_std; /* focal dist. (std) (cm) */ +extern float ct_PSpread_mean; /* pt. spread fn. (mean) (cm) */ +extern float ct_PSpread_std; /* pt. spread fn. (std) (cm) */ +extern float ct_Adjustment_std; /* adjustment dev. (std) (cm) */ +extern float ct_BlackSpot_rad; /* black spot radius (cm) */ +extern float ct_RMirror; /* rad. of single mirror (cm) */ +extern int ct_NMirrors; /* number of mirrors */ +extern float ct_CameraWidth; /* camera width (cm) */ +extern int ct_NPixels; /* number of pixels */ +extern float ct_PixelWidth; /* pixel width (cm) */ + +typedef struct +{ int i; /* mirror id */ + float f, /* focal length (cm) */ + sx, sy, /* curvilinear coords. (cm) */ + x, y, z, /* rectangular coords. (cm) */ + theta, phi, /* angles of mirror axis */ + xn, yn, zn; /* versor of mirror axis */ +} mirror; + +extern mirror *ct_data; /* ptr to mirror data */ +extern int nReflectivity; /* nr of refl. table elements */ +extern float *Reflectivity[]; /* ptr to refl. table */ +extern float *AxisDeviation[]; /* ptr to axisdev. table */ +extern float *ct_Focal; /* ptr to focal data */ + +/* Photon structures */ +typedef struct +{ float w, /* photon wavelength (nm) */ + x, y, /* (ground) impact point (cm) */ + u, v, /* direction cosines */ + t, /* arrival time (ns) */ + h; /* production height (cm) */ +} photon; + +typedef struct +{ float w, /* cphoton wavelength (nm) */ + x, y, /* (chamber) imp. point (cm) */ + u, v, /* direction cosines */ + t, /* arrival time (ns) */ + h, /* production height (cm) */ + phi; /* (chamber) inc. angle (rad) */ +} cphoton; +extern cphoton *CPhotons; + +/* Strings */ +#define RFL__START_MSG /* program, version */ \ + "\n ################################################\n" \ + " %s %s\n\n Simulation of the reflector\n" \ + " - J.C. Gonzalez - May 1998\n" \ + " Single file version\n" \ + " - D. Bastieri & C. Bigongiari - Jan 2000\n" \ + " ################################################\n\n" +#define SIGN_ERROR_FTL /* program, version */ \ + "Signature of parameters file is not correct.\n" \ + " *** Should be: **%s %s**.\n" +#define SKIP_TOKEN_MSG /* skipped token */ \ + " *** Skipping unknown token [%s].\n" +#define GEOM_NSPEC_FTL /* no parms */ \ + "No geometry file given.\n" \ + " *** Please specify its name after the \"ct_file\" directive.\n" +#define CPHS_ALLOC_FTL /* cphotons */ \ + "Cannot allocate enough memory space for %d cphotons.\n" +#define CPHS_ALLOC_LOG /* cphotons */ \ + "Allocated memory for %d cphotons.\n" +#define CERF_NFND__MSG /* cername */ \ + " *** Cerfile \"%s\" not found.\n" +#define CERF_OPEN__LOG /* cername */ \ + "Opened cerfile \"%s\".\n" +#define CERF_ERROR_LOG /* cername */ \ + "Error while reading cerfile \"%s\": skipping it.\n" +#define EVTN_WRONG_ERR /* first, last, cername */ \ + " *** Event boundaries [%ld, %ld] wrong for \"%s\".\n" \ + " *** Reading all events.\n" +#define TEMP_ERROR_FTL /* no parms */ \ + "Cannot open a temporary file: unable to handle overflow." +#define INFO_EVENT_LOG /* run, evt, energy */ \ + "Run %.0f, event %.0f: energy = %f.\n" +#define MEM__FREE__MSG /* no parms */ \ + "Freeing up allocated memory.\n" +#define VECT_FREE__LOG /* vector name */ \ + " ... freeing up [%s]\n" +#define RFLF_CLOSE_LOG /* no parms */ \ + "Closing reflector file.\n" +#define RFL__EXIT__MSG /* program, version */ \ + "%s %s: done.\n\n" + +#endif Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/input =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/input (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/input (revision 725) @@ -0,0 +1,31 @@ +reflector 0.4 +# +# Sample parameters file +# +verbose_level 1 +# +fixed_target 0. 90. +#range_events 1 2946 +#skip 1 2946 +#range_events 381 382 +max_events 1000000000 +#energy_cuts 800. 1000. +# +ct_file ../Data/magic.def +# +output_file my.rfl +# +atm_model ATM_CORSIKA +#atm_model ATM_NOATMOSPHERE +# +#data_paths 1 +#/scratch/work/cer000001 +# +#data_from_stdin +#data_to_stdout +# +#parallel_beam 0. 0. 0. 0. 2000. 2000. 300 300 1.e7 +#pm_parallel_beam unnamed.ppm 2. 10000. +# +cer_files +/big0/Maggi/MmcsData/cer014006 Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/lagrange.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/lagrange.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/lagrange.h (revision 725) @@ -0,0 +1,68 @@ +/* +///////////////////////////////////////////////////////////////// +// +// lagrange +//_______________________________________________________________ +// +// Created: Sun Jun 14 14:10:18 MET DST 1998 +// Author: Jose Carlos Gonzales +// Purpose: Macro for Lagrange interpolation of 3rd. order +// Notes: +// +///////////////////////////////////////////////////////////////// + + +//++ +// Formula for Lagrange interpolation of 3rd. order +// x: value to be interpolated +// t: table(2xN), table[0]: abscissas, table[1]: ordinates +// n: higher value of abscissas, such that t[0][n] <= x +//-- +*/ + +#define Lagrange(t,n,x) ((t[1][ (n) ]*((x-t[0][(n)+1])*(x-t[0][(n)+2]))/ \ + ((t[0][ (n) ]-t[0][(n)+1])*(t[0][ (n) ]-t[0][(n)+2])))+ \ + (t[1][(n)+1]*((x-t[0][ (n) ])*(x-t[0][(n)+2]))/ \ + ((t[0][(n)+1]-t[0][ (n) ])*(t[0][(n)+1]-t[0][(n)+2])))+ \ + (t[1][(n)+2]*((x-t[0][ (n) ])*(x-t[0][(n)+1]))/ \ + ((t[0][(n)+2]-t[0][ (n) ])*(t[0][(n)+2]-t[0][(n)+1]))) \ + ) + +/* +//++ +// Macro to find, and save in variable "m", the value of +// "n" to be used in the "Lagrange{t,n,x)" macro +//-- +*/ + +#define FindLagrange(t,m,x) {m = 0; while (t[0][++m] < x);} --m + +/* +////////////////////////////////////////////////////////////////////////////// +// Here follows a sample program using this two macros +////////////////////////////////////////////////////////////////////////////// +// #include +// #include "lagrange.h" +// +// void main(void) +// { +// float data[2][20]; +// int i, number; +// float x, y; +// +// for (i=0; i<20; ++i) { +// data[0][i] = i*10.; +// data[1][i] = 3.0*data[0][i]*data[0][i]; +// cout << data[0][i] << ' ' << data[1][i] << '\n'; +// } +// +// while (1==1) { +// cout << "Enter x = "; +// cin >> x; +// FindLagrange(data,number,x); +// y = Lagrange(data,number,x); +// cout << x << ' ' << y << '\n'; +// } +// } +////////////////////////////////////////////////////////////////////////////// */ + Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c (revision 725) @@ -0,0 +1,144 @@ +#include +#include +#include +#include "diag.h" +#include "parms.h" +#include "init.h" + +extern char whites[]; /* white chars (init) */ +extern char line[]; /* parsing buf. (init) */ + +/* Prototypes */ +extern void setall(long iseed1,long iseed2); /* rnds */ +static void ReadCerfiles(FILE *parfile); + +static void ReadCerfiles(FILE *parfile) +{ char *value_ptr = NULL; /* ptr at parm value */ + extern FILE *GetNextFile(char *cername); /* in main.c */ + + filelist = parfile; + parfile = NULL; + if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL || + (value_ptr=strtok(line, whites)) == NULL) + FatalError(FLST_NSPEC_FTL); + else if (value_ptr[0] == '@') + { fclose(filelist); + if ((filelist=fopen(value_ptr+1, "r")) == NULL) + FatalError(FLST_NFND__FTL, value_ptr+1); + else if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) + FatalError(FLST_NSPEC_FTL); + value_ptr = strtok(line, whites); } + + /* Set cername and find out event number bounds */ + strcpy(cername, value_ptr); + if ((value_ptr=strtok(NULL, whites)) == NULL) + { first_Event = 0; + last_Event = 1000000; } + else + { first_Event = atol(value_ptr); + value_ptr = strtok(NULL, whites); + last_Event = value_ptr ? atol(value_ptr) : 1000000; } + + /* Try to open cerfile. */ + + if ((cerfile=fopen(cername, "r")) == NULL) + { Message(CERF_NFND__MSG, cername); + cerfile=GetNextFile(cername); } + + /* If no valid cerfile is found then exit */ + if (cerfile == NULL) + FatalError(CERF_NSPEC_FTL); + + /* Check boundaries */ + if (first_Event > last_Event) + { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername); + first_Event = 0; + last_Event = 1000000; } + +} /* end of ReadCerfiles */ + +void ParmsSwitch(FILE *parfile) +{ char *value_ptr = NULL; /* ptr at parm value */ + int switch_end = FALSE; /* bool to exit loop */ + extern FILE *geofile; /* geo file (init) */ + extern void SetVerbose(int vlevel); /* from diag.c */ + extern void SetAtmModel(char *model); /* from atm.c */ + extern int ParseLine(FILE *parfile, /* FILE with parms */ + const char *token_list[], /* array w/tokens */ + int tokens, /* nr of tokens */ + char **value_ptr); /* ptr->parm val. */ + + do + { switch(ParseLine(parfile, parms, ARRAY_SZ(parms), &value_ptr)) + { case output_file: + if ((rflfile=fopen(value_ptr, "w+")) == NULL) + FatalError(OUTF_ERROR_FTL, value_ptr); + Message(OUTF_OPEN__MSG, value_ptr); + break; + case ct_file: + if ((geofile=fopen(value_ptr, "r")) == NULL) + FatalError(GEOF_ERROR_FTL, value_ptr); + Message(GEOF_OPEN__MSG, value_ptr); + strcat(strcpy(ct_BinaryName, value_ptr), ".mirr"); + ct_BinaryData = fopen(ct_BinaryName, "r"); + break; + case atm_model: + SetAtmModel(value_ptr); + break; + case verbose_level: + SetVerbose(atoi(value_ptr)); + break; + case fixed_target: + is_Fixed_Target = TRUE; + fixed_Theta = (float) atof(value_ptr); + value_ptr = strtok(NULL, whites); + if (value_ptr == NULL) + { Error(FIXD_TARGT_ERR); + is_Fixed_Target = FALSE; } + else + { fixed_Phi = (float) atof(value_ptr); + Message(FIXD_ENABL_MSG, fixed_Theta, fixed_Phi); + fixed_Theta *= (float) (M_PI/180.); + fixed_Phi *= (float) (M_PI/180.); } + break; + case max_events: + Message(MAX__EVTS__MSG, max_Events=atol(value_ptr)); + break; + case energy_cuts: + low_Ecut = (float) atof(value_ptr); + value_ptr = strtok(NULL, whites); + if (value_ptr == NULL) + { Error(ENRG_LIMIT_ERR); + low_Ecut = 0.; } + else + { high_Ecut = (float) atof(value_ptr); + Message(ENRG_CUTS__MSG, low_Ecut, high_Ecut); } + break; + case seeds: + Seeds[0] = atol(value_ptr); + value_ptr = strtok(NULL, whites); + if (value_ptr) Seeds[1] = atol(value_ptr); + else + { Error(SEED_ERROR_ERR); + Seeds[0] = 3141592L; } + break; + case random_pointing: + case repeat_random: +/********************************/ + break; + + case cer_files: + ReadCerfiles(parfile); + switch_end = TRUE; + default: switch_end = TRUE; + break; }} + while (!switch_end); + + if (filelist == NULL) + FatalError(FLST_NSPEC_FTL); + + /* Set random seeds */ + setall(Seeds[0], Seeds[1]); + Message(SEED_SET___MSG, Seeds[0], Seeds[1]); + +} /* end of ParmsSwitch */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h (revision 725) @@ -0,0 +1,61 @@ +#ifndef __RFL_PARMS__ +#define __RFL_PARMS__ + +/* token definition */ +#define ITEM_LIST /* LIST OF ITEMS IN THE PARAMETERS FILE */ \ +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), /* position towards which CT is pointing */ \ +T(max_events), /* maximum number of event to read */ \ +T(energy_cuts), /* lowest/highest energy allowed */ \ +T(seeds), /* seeds for random number generation */ \ +T(random_pointing), /* random CT pointing from each shower (hadrons) */ \ +T(repeat_random), /* number of times a random pointing is to be done */ \ +T(cer_files) /* start of filename list (must be last) */ + +#define T(x) x /* define T() as the name as it is */ + enum { ITEM_LIST }; +#undef T + +#define T(x) #x /* define T() as the string of x */ + const char *parms[] = { ITEM_LIST }; +#undef T +#undef ITEM_LIST + +/* Strings */ +#define OUTF_ERROR_FTL /* output filename */ \ + "Can't open output file: %s.\n" +#define OUTF_OPEN__MSG /* output filename */ \ + "Opened file \"%s\" as reflector file.\n" +#define GEOF_ERROR_FTL /* geometry filename */ \ + "Can't open geometry file: %s.\n" +#define GEOF_OPEN__MSG /* geometry filename */ \ + "Opened file \"%s\" as geometry file.\n" +#define FIXD_TARGT_ERR /* no parms */ \ + "Error while parsing \"fixed_target\" fields.\n" \ + " *** Reverted back to no \"fixed_target\" condition.\n" +#define FIXD_ENABL_MSG /* theta, phi */ \ + "Using \"fixed_target\" mode with theta=%.2fdeg and phi=%.2fdeg.\n" +#define MAX__EVTS__MSG /* max. nr. of evts. */ \ + "Processing at most %ld events.\n" +#define ENRG_LIMIT_ERR /* no parms */ \ + "Error while parsing \"energy_cuts\" fields.\n" \ + " *** Reverted back to no \"energy_cuts\" condition.\n" +#define ENRG_CUTS__MSG /* Elow, Ehigh */ \ + "Using \"energy_cuts\" mode with Elow=%.2fGeV and Ehigh=%.2fGeV.\n" +#define SEED_ERROR_ERR /* no parms */ \ + "Error while parsing second seed.\n" +#define SEED_SET___MSG /* Seed1, Seed2 */ \ + "Using seeds %ld and %ld.\n" +#define FLST_NSPEC_FTL /* no parms */ \ + "No input file specified.\n" \ + " *** Please specify some input file after the cer_files directive.\n" +#define FLST_NFND__FTL /* filelist */ \ + "Cannot find filelist \"%s\".\n" +#define CERF_NSPEC_FTL /* no parms */ \ + "No specified input file could be opened.\n" \ + " *** Please specify some *valid* filename after the cer_files directive.\n" + +#endif Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c (revision 725) @@ -0,0 +1,735 @@ +#include +#include +#include "diag.h" +#include "init.h" +#include "lagrange.h" +/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ + +/* random numbers */ +#define RandomNumber ranf() + +/* Speed of Light in vacuum, in m/s */ +#define Speed_of_Light_vacuum 299792458.0f +#define Speed_of_Light_air (Speed_of_Light_vacuum / 1.000293f) + +/* Speed of Light in vacuum, in cm/ns */ +#define Speed_of_Light_vacuum_cmns (Speed_of_Light_vacuum / 1.0e7f) +#define Speed_of_Light_air_cmns (Speed_of_Light_air / 1.0e7f) + +/* Macros */ +#define SQR(A) ((A)*(A)) +#define NORM(A) ((float) sqrt((SQR(A[0]))+(SQR(A[1]))+(SQR(A[2])))) + +/* Function declarations */ +extern float ranf(void); +void rnormal(double *r, int n); + +/* Static definitions */ +float OmegaCT[3][3]; +float OmegaICT[3][3]; +float Omega[3][3]; +float OmegaI[3][3]; + +static double NormalRandomNumbers[500]; +/* + From photons on ground, i.e. observation level, + to photons on focal plane, i.e. chamber ! + + Mirror reflectivity + Mirror reflection + Photon position on chamber + Position smearing + Timing + + Returned values: + + 0 OK photon reached the chamber + + 1 Photon lost due to mirror reflectivity + 2 Photon lost because out of mirror + 3 Photon lost due to black spot + 4 Photon lost because reflected out of chamber + +*/ + +int ph2cph(photon *ph, cphoton *cph) +{ + + float u, v, w; /* photon director cosines */ + + float r[3]; /* photon trajectory */ + float x[3]; /* position of the photon on ground */ + + float rCT[3]; /* photon trajectory in the system of the CT */ + float xCT[3]; /* photon position on ground (CT) */ + float rm[3]; /* photon trajectory in the system of a mirror */ + float xmm[3]; /* intermediate values */ + float xm[3]; /* photon position on ground */ + float xcut[3]; /* location of the cut sphere-trajectory */ + float xcutCT[3]; /* location of the cut sphere-trajectory (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 trajectory-mirror center */ + + float wl; /* photon wavelength */ + float reflec; /* reflectivity for a photon */ + + float h; /* photon production height */ + + int i, k; /* simple counters */ + int i_mirror=-1; /* number of a given mirror */ + + + float distmirr, distmirr2; /* distances used in MAGIC reflection routine */ + float sx, sy; + float t1, t2; + + void makeOmega(float theta, float phi); + void makeOmegaI(float theta, float phi); + void applyMxV(float M[3][3], float *V, float *Vp); + float Lin2Curv(float x); + + /* begin code */ + + /* get photon wawelength */ + + wl = ph->w; + + /* get position on ground */ + + x[0] = ph->x; + x[1] = ph->y; + x[2] = 0.0; /* ground => obs. level => z=0 */ + + /* get director cosines x,y on ground */ + + r[0] = ph->u; + r[1] = ph->v; + r[2] = (float) sqrt(1.0 - r[0]*r[0] - r[1]*r[1]); + + /* get photon time and production height */ + + h = ph->h; + + /*!@' + + @#### 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: 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 ) return 1; + + + /*!@' + + @#### Reflection on mirrors. + We calculate reflected photon direction + + */ + + /* ++ + REFLECTION + -- */ + + Debug("@1 %f %f %f %f %f %f\n", x[0], x[1], x[2], r[0], r[1], r[2]); + + /* change to the system of the CT */ + + applyMxV( OmegaCT, x, xCT ); + applyMxV( OmegaCT, r, rCT ); + + /* + before moving to the system of the mirror + we look whether the photon hits a mirror or not + + calculate the intersection of the trajectory 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 = (float) (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 = (float) sqrt( b*b - 4.0*a*c ); + + /* two possible values for z */ + + t1 = (float) ((-b+d) / (2.0*a)); + t2 = (float) ((-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 trajectory 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]); + + /* 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 << "CONDITION 1 !" << endl << flush; + cout << '1'; + */ + return 2; + } + + /* calculate the mirror to be used */ + + distmirr = 1000000.0f; + + for (i=0; i=ct_RMirror; ++i) { + distmirr2 = (float) sqrt(SQR(ct_data[i].x - xcut[0]) + + SQR(ct_data[i].y - xcut[1]) + + SQR(ct_data[i].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 ((fabs(ct_data[i_mirror].sx - sx) > ct_RMirror) || + (fabs(ct_data[i_mirror].sy - sy) > ct_RMirror)) { + /* + cout << "CONDITION 2 !" << endl << flush; + cout << '2'; + */ + return 2; + } + + /* calculate matrices for the mirror */ + + makeOmega (-ct_data[i_mirror].theta, + ct_data[i_mirror].phi); + makeOmegaI(-ct_data[i_mirror].theta, + ct_data[i_mirror].phi); + + /* change to the system of the mirror */ + + /* first translation... */ + + xmm[0] = xCT[0] - ct_data[i_mirror].x; + xmm[1] = xCT[1] - ct_data[i_mirror].y; + xmm[2] = xCT[2] - ct_data[i_mirror].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 trajectory 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 coefficient + shown above + */ + + a = SQR(rm[0]) + SQR(rm[1]) + SQR(rm[2]); + + b = (float) (-2*(2.*ct_data[i_mirror].f*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 = (float) sqrt( b*b - 4.0*a*c ); + + /* two possible values for z */ + + t1 = (float) ((-b+d) / (2.0*a)); + t2 = (float) ((-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'; + */ + return 3; + } + + /* + if we still have the photon, we continue with the reflexion; + we calculate normal vector in this point + (and normalize, with the sign changed) + */ + + rnor[0] = 2.0f*xcut[0]; + rnor[1] = 2.0f*xcut[1]; + rnor[2] = (float) (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 trajectory 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 = (float) fabs(rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]); + + /* finally!!! we have the reflected trajectory of the photon */ + + rrefl[0] = (float) (2.0*rnor[0]*calpha - rm[0]); + rrefl[1] = (float) (2.0*rnor[1]*calpha - rm[1]); + rrefl[2] = (float) (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].x; + xcutCT[1] += ct_data[i_mirror].y; + xcutCT[2] += ct_data[i_mirror].z; + + /* + calculate intersection of this trajectory 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] += (float) (NormalRandomNumbers[0] * ct_PSpread_mean); + xcam[1] += (float) (NormalRandomNumbers[1] * ct_PSpread_mean); + + /* check whether the photon goes out of the camera */ + + if ( (SQR(xcam[0])+SQR(xcam[1])) > SQR(ct_CameraWidth) ) { + return 4; + } + + /* + ++ + 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 = (float) asin(rreflCT[2]); + + /* + ++ + TIMING + -- + */ + + /* calculate the new time of the photon (in the camera) */ + + t = ph->t; + + /* + substract path from the mirror till the ground, 'cos + the photon actually hit the mirror!! + */ + + t = (float) (t + ((( xm[2] > 0. ) ? -1.0 : +1.0) * + 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 = (float) (t + sqrt( SQR(xcutCT[0] - xcam[0]) + + SQR(xcutCT[1] - xcam[1]) + + SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns); + + /* show it */ + + Debug("@2 %f %f %f %f %f %f %f %f %f\n" + "@3 %f %f %d %f %f %f %f\n" + "@4 %f %f %f %f\n", + xCT[0], xCT[1], xCT[2], rCT[0], rCT[1], rCT[2], + xcut[0], xcut[1], xcut[2], + sx, sy, i_mirror, ct_data[i_mirror].sx, ct_data[i_mirror].sy, + ct_data[i_mirror].sx - sx, ct_data[i_mirror].sy - sy, + xcam[0], xcam[1], xcam[2], phi); + + /* Output */ + +/* cph->w = wl; */ + cph->x = xcam[0]; + cph->y = xcam[1]; + cph->u = r[0]; + cph->v = r[1]; + cph->t = t; + cph->h = h; + cph->phi = phi; + + return 0; + +} /* end of ph2cph */ + + +/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +!--------------------------------------------------------------------- + @name makeOmega + + @desc function to calculate the matrix Omega(theta,phi) + + @var theta Angle theta of the transformation + @var phi Angle phi of the transformation + + @date Sat Jun 27 05:58:56 MET DST 1998 +---------------------------------------------------------------------- + @function +*/ + +void +makeOmega (float theta, float phi) +{ + static float ct, st, cp, sp; + + /* shortcuts for cosine and sine of theta and phi */ + ct = (float) cos(theta); + st = (float) sin(theta); + cp = (float) cos(phi); + sp = (float) sin(phi); + + /* save values in the array (see top of file) */ + Omega[0][0] = cp*ct; + Omega[0][1] = sp*ct; + Omega[0][2] = -st; + + Omega[1][0] = -sp; + Omega[1][1] = cp; + Omega[1][2] = 0; + + Omega[2][0] = cp*st; + Omega[2][1] = sp*st; + Omega[2][2] = ct; +} + + +/* +!--------------------------------------------------------------------- + @name makeOmegaI + + @desc function to calculate the matrix Omega-1(theta,phi) + + @var theta Angle theta of the transformation + @var phi Angle phi of the transformation + + @date Sat Jun 27 05:58:56 MET DST 1998 +---------------------------------------------------------------------- + @function +*/ + +void +makeOmegaI(float theta, float phi) +{ + static float ct, st, cp, sp; + + /* shortcuts for cosine and sine of theta and phi */ + ct = (float) cos(theta); + st = (float) sin(theta); + cp = (float) cos(phi); + sp = (float) sin(phi); + + /* save values in the array (see top of file) */ + OmegaI[0][0] = cp*ct; + OmegaI[0][1] = -sp; + OmegaI[0][2] = cp*st; + + OmegaI[1][0] = sp*ct; + OmegaI[1][1] = cp; + OmegaI[1][2] = sp*st; + + OmegaI[2][0] = -st; + OmegaI[2][1] = 0; + OmegaI[2][2] = ct; +} + + +/* +!--------------------------------------------------------------------- + @name applyMxv + + @desc returns the vector v' such that v' = M x v + + @var M matrix of the transformation + @var v vector to be multiplied + @var vi resulting vector + + @date Sat Jun 27 05:58:56 MET DST 1998 +---------------------------------------------------------------------- + @function +*/ + +void +applyMxV(float M[3][3], float *V, float *Vp) +{ + Vp[0] = (M[0][0] * V[0] + + M[0][1] * V[1] + + M[0][2] * V[2]); + Vp[1] = (M[1][0] * V[0] + + M[1][1] * V[1] + + M[1][2] * V[2]); + Vp[2] = (M[2][0] * V[0] + + M[2][1] * V[1] + + M[2][2] * V[2]); +} + +/* +!--------------------------------------------------------------------- + @name Lin2Curv + + @desc Linear (Euclidean) to Curvilinear distance + + @var x Radial distance from the axis of the paraboloid + + @return Curvilinear distance over the parabolic shape + + @date Wed Jul 8 15:25:39 MET DST 1998 +---------------------------------------------------------------------- + @function +*/ + +float +Lin2Curv(float x) +{ + x /= 100.f; + return ((x + 0.000144175317185f * x * x * x)*100.f); +} + +/*!--------------------------------------------------------------------- +// @name rnormal +// +// @desc returns n(=2k) normaly distributed numbers +// +// @var *r pointer to a vector where we write the numbers +// @var n how many numbers do we generate +// +// @date Sat Jun 27 05:58:56 MET DST 1998 +//---------------------------------------------------------------------- +// @function */ + +void rnormal(double *r, int n) +{ + + double z1, z2; + int i; + + for (i=0; i +#include +#include +#include +#include "version.h" +#include "diag.h" +#include "init.h" +#include "header.h" + +FILE *rflfile = NULL, /* reflector (output) file */ + *cerfile = NULL; /* current cerfile */ + +cphoton *CPhotons = NULL; /* cphoton array */ +static photon Photons[PH_IN_DATABLOCK]; /* photon array */ + +static long readp = 0L; /* read ptr (on cerfile) */ +static long writep = 0L; /* write ptr (on rflfile) */ + + +extern CerHeader *cheadp; /* var inited in header.c */ +extern RflHeader *rheadp; /* var inited in header.c */ + +/* Prototypes */ +void clean(void); +FILE *GetNextFile(char *cername); +static int GetEvent(void); +static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile); + +void main(void) +{ long event = 0L; /* event counter */ + + /* Read init & geometry parms, init files and vars */ + init(NULL); + + /* Processing loop */ + while(event < max_Events && GetEvent()) + { if (ProcessEvent(cheadp, cerfile, rflfile)) event++; + if (event % SHOW_ME == 0) Log(INFO_EVENT_LOG, + cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal); } + + /* Writing final flags */ + fwrite(FLAG_END_OF_FILE, SIZE_OF_FLAGS, 1, rflfile); + + /* Close file */ + Log(RFLF_CLOSE_LOG); + fclose(rflfile); + + /* Clean memory and exit */ + clean(); + Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION)); + +} /* end of main */ + +static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile) +{ extern int absorption(float wlen, float height, float theta); + extern int ph2cph(photon *ph, cphoton *cph); + extern float Omega[3][3]; + extern float OmegaI[3][3]; + extern float OmegaCT[3][3]; + extern float OmegaICT[3][3]; + extern void makeOmega(float theta, float phi); + extern void makeOmegaI(float theta, float phi); + + /* Various counters: phs = absphs + refphs[0..3] + cphs */ + long phs, /* Number of incoming photons */ + absphs, /* Photons absorbed */ + refphs[4], /* Photons not reflected */ + cphs; /* Number of cphotons */ + + FILE *tmpf=NULL; /* Temp fp to handle o/f */ + size_t read; /* items read: != 1 on error */ + int overflow=0, /* how many o/f on cphs */ + ref_type, /* ret value from reflection */ + ph; /* photon number */ + + float first = 1e8f, /* Photon arrival times */ + last = 0; + + /* photon quantities */ + float wlen, /* photon wavelength */ + theta; /* photon zenith angle */ + + /* Reset counters and set fileptr */ + phs= absphs= refphs[0]= refphs[1]= refphs[2]= refphs[3]= cphs= 0; + writep = ftell(rflfile); + + /* Translate the EVTH from cerfile to rflfile */ + TranslateHeader(rheadp, cheadp); + + /* Calculate OmegaCT matrices */ + if (is_Fixed_Target) + { makeOmega(fixed_Theta, fixed_Phi); + makeOmegaI(fixed_Theta, fixed_Phi); } + else + { makeOmega(cheadp->Theta, cheadp->Phi); + makeOmegaI(cheadp->Theta, cheadp->Phi); } + memcpy( OmegaCT, Omega, 9*sizeof(float) ); + memcpy( OmegaICT, OmegaI, 9*sizeof(float) ); + + /* Loop on data blocks */ + while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile)) + && strncmp((char *)Photons, "EVTE", 4)) + + /* Loop on phs inside block: exit when wlength is 0 */ + { for (ph=0; phCorePos[0][0]; + Photons[ph].y -= cheadp->CorePos[1][0]; + + /* Increment number of photons */ + phs++; + + /* Calculate some quantities */ + theta = (float) acos(sqrt( + 1.f - Photons[ph].u*Photons[ph].u + - Photons[ph].v*Photons[ph].v)); + /* Check absorption */ + + if (absorption(wlen, Photons[ph].h, theta)) + absphs++; + + /* Check reflection */ + else if (0 != (ref_type = + ph2cph(&Photons[ph], &CPhotons[cphs]))) + refphs[ref_type-1]++; + + else /* Photon passed */ + { + /* Update first/last arrival times */ + if (first > CPhotons[cphs].t) first = CPhotons[cphs].t; + if ( last < CPhotons[cphs].t) last = CPhotons[cphs].t; + + /* Update cphs */ + if (++cphs == NR_OF_CPHOTONS) /* Overflow */ + { + /* if it is the first o/f open a tempfile */ + if (overflow++ == 0) + { Log("Spooling... "); + tmpf = tmpfile(); + if (tmpf == NULL) FatalError(TEMP_ERROR_FTL); } + + /* Write now the cphoton array */ + fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf); + + /* Reset the cphs counter */ + cphs = 0; + + } /* if overflow */ + } /* else (=Photon passed) */ + } /* end of loop inside datablock */ + } /* end of loop on datablocks */ + + /* Check if there was an error while reading cerfile */ + if (read != 1) fseek(rflfile, writep, SEEK_SET); + + else /* no error: write the new event */ + + { /* Write "start of event" flag */ + fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile); + + /* Update arrival times */ + rheadp->TimeFirst = first; + rheadp->TimeLast = last; + + /* Update RflHeader with info on ph/cph nrs and write it */ + rheadp->CORSIKAPhs = phs; + rheadp->AtmAbsPhs = absphs; + rheadp->MirrAbsPhs = refphs[0]; + rheadp->OutOfMirrPhs = refphs[1]; + rheadp->BlackSpotPhs = refphs[2]; + rheadp->OutOfChamPhs = refphs[3]; + rheadp->CPhotons = (long) overflow * NR_OF_CPHOTONS + cphs; + fwrite(rheadp, sizeof(RflHeader), 1, rflfile); + + /* If there was an overflow, append data from tempfile */ + if (overflow) + { + /* Unload data from CPhotons */ + fwrite(CPhotons, sizeof(cphoton), cphs, tmpf); + + /* Transfer data */ + fseek(tmpf, 0L, SEEK_SET); /* Start from the beginning */ + while (overflow--) + { fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf); + fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile); } + + /* Reload data in CPhotons */ + fread (CPhotons, sizeof(cphoton), cphs, tmpf); + + /* Close (and delete) temp file */ + fclose(tmpf); } + + /* Write (remaining) cphotons */ + fwrite(CPhotons, sizeof(cphoton), cphs, rflfile); + + /* Write "end of event" flag */ + fwrite(FLAG_END_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile); + + Debug(">>> %6ld = %6ld + %6ld + %6ld + %6ld + %6ld + %6f\n", + phs, absphs, + refphs[0], refphs[1], refphs[2], refphs[3], + rheadp->CPhotons); } + + return read == 1; +} /* end of ProcessEvent */ + +static int GetEvent(void) +{ int found = FALSE, /* event found */ + isWrong = FALSE; /* cerfile is wrong */ + static int newFile = TRUE; /* if TRUE, check if cerfile is valid */ + + do + { /* In case the just-opened file is a valid cerfile, + starting with a RUNH, loop until a valid event is found: + id est with: first_Event <= EvtNumber <= last_Event + and low_Ecut <= Etotal <= high_Ecut + If the search was successful, "found" is set to TRUE. + If there are reading errors, "isWrong" is set to TRUE. */ + + if (newFile + && (1 != fread(cheadp, sizeof(CerHeader), 1, cerfile) + || strncmp(cheadp->EVTH, "RUNH", 4))) + { isWrong = TRUE; } + + else do + { /* Push fileptr in case it is a "EVTH" */ + readp = ftell(cerfile); + + /* Error: exit loop */ + if (1 != fread(cheadp, sizeof(CerHeader), 1, cerfile)) + { isWrong = TRUE; break; } + + /* Ok: set found at TRUE and exit loop */ + if (strncmp(cheadp->EVTH, "EVTH", 4) == 0 + && first_Event <= cheadp->EvtNumber + && cheadp->EvtNumber <= last_Event + && low_Ecut <= cheadp->Etotal + && cheadp->Etotal <= high_Ecut) + { found = TRUE; } + + /* File is finished: exit loop */ + else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0) + break; + + } while(found == FALSE); + + /* If there was an error, log it */ + if (isWrong) Log(CERF_ERROR_LOG, cername); + + /* If a new event was not found, get, if any, a new cerfile. + "newFile" is set in this case, to check the validity of + just-opened cerfile. If GetNextFile return a NULL (no + new cerfile), then exit thoroughly. */ + + if (found) newFile = FALSE; + else + { if ((cerfile=GetNextFile(cername)) == NULL) break; + newFile = TRUE; } + + } while(found == FALSE); + + return found; +} /* end of GetEvent */ + +/* Files are all CER-like files (cerfiles). + Filenames are written one per line and may be followed + by a ",FIRST_EVT,LAST_EVT" directive, such as: + -> cer19999 15 167 + telling to deal with events from 15 to 167 in file "cer19999" + The file list may be appended to the input file or can be stored + in a separate file pointed by an '@' directive contained in the + input file. In both cases, file referral follows the tag + "cer_files". Files not accessible are thoroughly skipped. + GetNextFile gets and opens the next readable cerfile. */ + +FILE *GetNextFile(char *cername) +{ FILE *f = NULL; /* return value (cerfile ptr) */ + char *value_ptr; /* ptr at parm value */ + extern char line[]; /* white chars (init) */ + extern char whites[]; /* parsing buf. (init) */ + + /* Close, if any, the previous cerfile, writing "END_OF_RUN". */ + if (cerfile) + { fclose(cerfile); + fwrite(FLAG_END_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); } + + /* Parse next line in filelist */ + do + { if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) break; + if ((value_ptr=strtok(line, whites)) == NULL) continue; + + /* If you found a line with some meaning, try to open the file */ + if ((f=fopen(value_ptr, "r")) == NULL) + Message(CERF_NFND__MSG, value_ptr); + + else /* Cerfile is readable */ + { /* Store its name in "cername" */ + Log(CERF_OPEN__LOG, strcpy(cername, value_ptr)); + + /* Read, if any, event bounds */ + if ((value_ptr=strtok(NULL, whites)) == NULL) + { first_Event = 0; + last_Event = 1000000; } + else + { first_Event = atol(value_ptr); + value_ptr = strtok(NULL, whites); + last_Event = value_ptr ? atol(value_ptr) : 1000000; } + + /* Check boundaries */ + if (first_Event > last_Event) + { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername); + first_Event = 0; + last_Event = 1000000; }}} + + while (f == NULL); /* Loop until a readable file is found */ + + /* If a new cerfile is met, write "START_OF_RUN" */ + if (f) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); + + return f; +} /* end of GetNextFile */ + +/* This function frees up allocated memory before exiting */ + +void clean(void) +{ Message(MEM__FREE__MSG); + + if (ct_data) + { free(ct_data); Log(VECT_FREE__LOG, "ct_data"); } + + if (Reflectivity[0]) + { free(Reflectivity[0]); Log(VECT_FREE__LOG, "Reflectivity[0]"); } + if (Reflectivity[1]) + { free(Reflectivity[1]); Log(VECT_FREE__LOG, "Reflectivity[1]"); } + + if (AxisDeviation[0]) + { free(AxisDeviation[0]); Log(VECT_FREE__LOG, "AxisDeviation[0]"); } + if (AxisDeviation[1]) + { free(AxisDeviation[1]); Log(VECT_FREE__LOG, "AxisDeviation[1]"); } + + if (ct_Focal) + { free(ct_Focal); Log(VECT_FREE__LOG, "ct_Focal"); } + + if (CPhotons) + { free(CPhotons); Log(VECT_FREE__LOG, "CPhotons"); } +} /* end of clean */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/version.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/version.h (revision 725) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/version.h (revision 725) @@ -0,0 +1,10 @@ +#ifndef __RFL_VERSION__ +#define __RFL_VERSION__ + +#define QUOTE_TP(x) #x +#define QUOTE(x) QUOTE_TP(x) + +#define PROGRAM reflector +#define VERSION 0.4 + +#endif