Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile (revision 1431) @@ -19,7 +19,7 @@ # # $RCSfile: Makefile,v $ -# $Revision: 1.2 $ -# $Author: blanch $ -# $Date: 2001-04-26 07:57:41 $ +# $Revision: 1.3 $ +# $Author: bigongia $ +# $Date: 2002-07-24 14:35:45 $ # ################################################################## @@ -53,9 +53,8 @@ # compilation and linking flags -CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} ${DEBUG} +#CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} ${DEBUG} +CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} CFLAGS = ${CXXFLAGS} -#CFLAGS = -std1 -warnprotos -msg_enable level6 -msg_disable returnchecks FFLAGS = ${CXXFLAGS} -LIBS = ${RANLIB} ${FORLIBS} LIBS = ${CERNLIB} ${RANLIB} ${FORLIBS} @@ -180,9 +179,5 @@ diag.o: /usr/include/bits/types.h /usr/include/libio.h diag.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h -diag.o: /usr/include/stdlib.h /usr/include/sys/types.h /usr/include/time.h -diag.o: /usr/include/endian.h /usr/include/bits/endian.h -diag.o: /usr/include/sys/select.h /usr/include/bits/select.h -diag.o: /usr/include/bits/sigset.h /usr/include/sys/sysmacros.h -diag.o: /usr/include/alloca.h version.h diag.h +diag.o: /usr/include/stdlib.h version.h diag.h init.o: /usr/include/stdio.h /usr/include/features.h /usr/include/sys/cdefs.h init.o: /usr/include/gnu/stubs.h @@ -191,13 +186,7 @@ init.o: /usr/include/bits/types.h /usr/include/libio.h init.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h -init.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/sys/types.h -init.o: /usr/include/time.h /usr/include/endian.h /usr/include/bits/endian.h -init.o: /usr/include/sys/select.h /usr/include/bits/select.h -init.o: /usr/include/bits/sigset.h /usr/include/sys/sysmacros.h -init.o: /usr/include/alloca.h /usr/include/math.h +init.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/math.h init.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h -init.o: /usr/include/bits/mathcalls.h -init.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h -init.o: version.h diag.h init.h +init.o: /usr/include/bits/mathcalls.h version.h diag.h init.h parms.o: /usr/include/stdio.h /usr/include/features.h parms.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h @@ -206,9 +195,7 @@ parms.o: /usr/include/bits/types.h /usr/include/libio.h parms.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h -parms.o: /usr/include/string.h /usr/include/math.h +parms.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/math.h parms.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h -parms.o: /usr/include/bits/mathcalls.h -parms.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h -parms.o: diag.h parms.h init.h +parms.o: /usr/include/bits/mathcalls.h diag.h parms.h init.h geometry.o: /usr/include/stdio.h /usr/include/features.h geometry.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h @@ -217,9 +204,7 @@ geometry.o: /usr/include/bits/types.h /usr/include/libio.h geometry.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h -geometry.o: /usr/include/string.h /usr/include/math.h +geometry.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/math.h geometry.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h -geometry.o: /usr/include/bits/mathcalls.h -geometry.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h -geometry.o: diag.h geometry.h init.h +geometry.o: /usr/include/bits/mathcalls.h diag.h geometry.h init.h atm.o: /usr/include/stdio.h /usr/include/features.h /usr/include/sys/cdefs.h atm.o: /usr/include/gnu/stubs.h @@ -229,7 +214,6 @@ atm.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h atm.o: /usr/include/string.h /usr/include/math.h /usr/include/bits/huge_val.h -atm.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h -atm.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h diag.h -atm.o: atm.h init.h +atm.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h diag.h atm.h +atm.o: init.h ph2cph.o: /usr/include/stdio.h /usr/include/features.h ph2cph.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h @@ -239,7 +223,6 @@ ph2cph.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h ph2cph.o: /usr/include/math.h /usr/include/bits/huge_val.h -ph2cph.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h -ph2cph.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h -ph2cph.o: diag.h init.h lagrange.h +ph2cph.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h diag.h +ph2cph.o: init.h lagrange.h header.o: /usr/include/string.h /usr/include/features.h header.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h @@ -252,12 +235,5 @@ reflector.o: /usr/include/bits/types.h /usr/include/libio.h reflector.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h -reflector.o: /usr/include/stdlib.h /usr/include/sys/types.h -reflector.o: /usr/include/time.h /usr/include/endian.h -reflector.o: /usr/include/bits/endian.h /usr/include/sys/select.h -reflector.o: /usr/include/bits/select.h /usr/include/bits/sigset.h -reflector.o: /usr/include/sys/sysmacros.h /usr/include/alloca.h -reflector.o: /usr/include/string.h /usr/include/math.h +reflector.o: /usr/include/stdlib.h /usr/include/string.h /usr/include/math.h reflector.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h -reflector.o: /usr/include/bits/mathcalls.h -reflector.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h -reflector.o: version.h diag.h init.h header.h +reflector.o: /usr/include/bits/mathcalls.h version.h diag.h init.h header.h Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c (revision 1431) Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f (revision 1431) @@ -360,4 +360,13 @@ x4 = sqrt((2. * h + rt) / (2. * hscale)) +c-- AM Dec 2001, to avoid crash! A few photons seem to be "corrupted" +c-- (have absurd value) in a cer file... + + if (abs(x3-x4) .lt. 1.e-10) then + tr_atmos = -1. + RETURN + endif + + e1 = derfc(x1) e2 = derfc(x2) @@ -366,4 +375,5 @@ m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4)) + ********************************************************************** @@ -442,6 +452,7 @@ CON_OZ = 2 + 2001 IF (LONG(CON_OZ).GT.wavelength) THEN - A = (OZ_ABI(ROW,CON_OZ)-OZ_ABI(ROW,CON_OZ-1))/ + 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)) Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.c (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.c (revision 1431) @@ -8,10 +8,12 @@ extern char line[]; /* parsing buf. (init) */ +extern char axisdev_filename[256], reflectivity_filename[256]; + +float mean_refl; /* Mirror mean reflectivity 270-610 nm. + * AM June 2002. + */ 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) */ @@ -43,56 +45,64 @@ 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; }} + { + if (line[0] == '#') continue; + + if (nReflectivity == 0) + { + nReflectivity = atoi(line); + if (nReflectivity) + { + if ((Reflectivity[0] = + (float *) malloc(sizeof(float) * nReflectivity)) == NULL || + (Reflectivity[1] = + (float *) malloc(sizeof(float) * nReflectivity)) == NULL) + FatalError(REFL_ALLOC_FTL, nReflectivity); + } + } + else if (2 == sscanf(line, "%f %f", &Reflectivity[0][current], + &Reflectivity[1][current])) + { + // Added June 2002, AM: + mean_refl += Reflectivity[1][current]; + current++; + if (current >= nReflectivity) break; + } + } fclose(datfile); nReflectivity = current; + if (current > 0) + mean_refl /= (float) current; + } /* end of ReadReflectivity */ @@ -103,4 +113,7 @@ if (datfile == NULL) FatalError(AXIS_ERROR_FTL, datname); + else + printf("Reading file %s\n", axisdev_filename); + if ((AxisDeviation[0]= (float *) malloc(sizeof(float) * ct_NMirrors)) == NULL @@ -149,22 +162,16 @@ break; case focal_distance: - Log(LOG__FLOAT_LOG, "focal distance (average)", + Log(LOG__FLOAT_LOG, "focal distance (average, cm)", 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)); + case focal_std: /* not implemented. */ break; case point_spread: - Log(LOG__FLOAT_LOG, "point spread fn. (average)", + Log(LOG__FLOAT_LOG, "point spread fn. sigma (average, cm)", 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)); + case point_std: /* not implemented */ + break; + case adjustment_dev: /* not implemented */ break; case black_spot: @@ -172,14 +179,14 @@ ct_BlackSpot_rad = (float) atof(value_ptr)); break; + case n_mirrors: + Log(LOG__INT___LOG, "number of mirrors", + ct_NMirrors = atoi(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)", + Log(LOG__FLOAT_LOG, "camera radius (cm)", ct_CameraWidth = (float) atof(value_ptr)); break; @@ -191,4 +198,10 @@ Log(LOG__FLOAT_LOG, "pixel width (cm)", ct_PixelWidth = (float) atof(value_ptr)); + break; + case n_centralpixels: + /* this parameter is for camera, not for reflector */ + break; + case n_gappixels: + /* this parameter is for camera, not for reflector */ break; case refl_file: @@ -208,6 +221,14 @@ fclose(geofile); - if (Reflectivity[0] == NULL) ReadReflectivity(REFLECTIVITY_FILE); - if (AxisDeviation[0]== NULL) ReadAxisDev(AXISDEVIATION_FILE); + if (strlen(reflectivity_filename) == 0) + strcpy(reflectivity_filename, REFLECTIVITY_FILE); + if (Reflectivity[0] == NULL) ReadReflectivity(reflectivity_filename); + + if (strlen(axisdev_filename) == 0) + strcpy(axisdev_filename, AXISDEVIATION_FILE); + if (AxisDeviation[0]== NULL) ReadAxisDev(axisdev_filename); + + ReadFocals(); } /* end of GeometrySwitch */ + Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.h (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.h (revision 1431) @@ -17,6 +17,8 @@ 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 */ - +T(define_mirrors), /* this entry is followed by the def. of pixels */ \ +T(n_centralpixels), /* this token is not for reflector but for camera */ \ +T(n_gappixels) /* same comment as for previous token. */ + #define T(x) x /* define T() as the name as it is */ enum { ITEM_LIST }; @@ -48,12 +50,7 @@ "Located and opened file \"%s\" containing mirror data.\n" #define READ_ASCII_LOG /* no parms */ \ - "Start reading ASCII data for mirror nr:\n" + "Reading ASCII data for mirrors.\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" @@ -82,4 +79,6 @@ */ +extern float mean_refl; /* mean mirror reflectivity AM June 2002 */ + #define RANDOM_POINTING_MAX_SEPARATION 0.104719755119660 Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/header.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/header.c (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/header.c (revision 1431) @@ -9,35 +9,62 @@ static CerHeader chead; CerHeader *cheadp = &chead; +extern float fixed_Phi, fixed_Theta; +extern int ct_NMirrors; +extern float mean_refl; + 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; + 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->Corsika_version= c->Corsika_version; + r->NumObsLev = c->NumObsLev; + r->HeightLev = c->HeightLev[0]; - 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; + 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; -} /* end of TranslateHeader */ + memcpy(r->p, c->p, 3*sizeof(float)); + memcpy(r->RndData, c->RndData, 30*sizeof(float)); + memcpy(r->CorePos, c->CorePos, 40*sizeof(float)); + + /* Next 4 variables added in June 2002, AM */ + r->telescopePhi = fixed_Phi; + r->telescopeTheta = fixed_Theta; + r->num_mirrors = ct_NMirrors; + r->mean_reflectivity = mean_refl; + + if (c->Corsika_version >= 6.) /* Viewcone option implemented only in c6xx */ + { + if (c->viewcone_angles[0] > 0.) + { + printf ("ERROR: Input cer file was run using Corsika option\n"); + printf (" VIEWCONE VUECON(1) VUECON(2) with VUECON(1) > 0\n"); + printf (" The reflector program only supports VUECON(1) = 0\n\n"); + exit(-1); + } + else + r->ViewConeRadius = c->viewcone_angles[1]; /* degrees */ + } + else + r->ViewConeRadius = 0.; + +} /* end of TranslateHeader */ + + Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h (revision 1431) @@ -19,5 +19,5 @@ float RunNumber; float DateRun; - float VersionPGM; + float Corsika_version; float NumObsLev; @@ -72,8 +72,9 @@ float CorePos[2][20]; - float dmmy1; - float SpinTheta; - float SpinPhi; - float dmmy2[132]; + float dmmy1[14]; + float viewcone_angles[2]; /* (degrees) Inner and outer angles in + * Corsika's VIEWCONE option + */ + float dmmy2[119]; } CerHeader; @@ -94,8 +95,19 @@ float RunNumber; float DateRun; - float VersionPGM; + float Corsika_version; - float NumObsLev; - float HeightLev[10]; + float NumObsLev; /* Should be 1 for MAGIC simulation */ + float HeightLev; /* Observation Level */ + + /* Changed meaning of next 9 variables. June 2002, A.Moralejo: */ + float num_mirrors; + float mean_reflectivity; + float longi_Nmax; + float longi_t0; + float longi_tmax; + float longi_a; + float longi_b; + float longi_c; + float longi_chi2; float SlopeSpec; @@ -114,8 +126,22 @@ float TimeLast; - float deviationPhi; - float deviationTheta; + /* AM, 23/05/2002: Changed meaning of following + * three variables (which were unused before): + */ + float telescopePhi; /* rad */ + float telescopeTheta; /* rad */ - float Trigger; + float ViewConeRadius; /* Degrees. + * Radius of "view cone" when the primaries' + * directions generated by Corsika lie within + * a cone around a fixed direction. This is + * only possible with Corsika>6 versions. In + * that case, PhiMin=PhiMax and + * ThetaMin=ThetaMax (also in this header) + * indicate the axis of this cone. + * If ViewConeRadius==0, it means that + * the VIEWCONE option was not used. + */ + float CORSIKAPhs; /* Original photons written by Corsika */ @@ -124,7 +150,21 @@ 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 */ + float OutOfChamPhs; /* Photons outside the camera */ + float CPhotons; /* Photons reaching the camera */ } RflHeader; + + +typedef struct /* EVTE from cerfile */ +{ + float dummy1[255]; + float longi_Nmax; + float longi_t0; + float longi_tmax; + float longi_a; + float longi_b; + float longi_c; + float longi_chi2; + float dummy2[11]; +} Event_end; void TranslateHeader(RflHeader *r, CerHeader *c); Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c (revision 1431) @@ -20,4 +20,7 @@ float fixed_Theta, /* zenith angle (rad) */ fixed_Phi; /* azi (0N->90E) (rad) */ + +float Telescope_x = 0.; +float Telescope_y = 0.; /* Telescope coordinates (cm) */ int is_Random_Pointing=FALSE; /* random pointing? */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h (revision 1431) @@ -4,5 +4,5 @@ /* Constant definitions */ #define LINE_MAX_LENGTH 128 /* parsing buffer sz. */ -#define NR_OF_CPHOTONS 65535 /* CPhotons 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 */ @@ -41,4 +41,6 @@ extern float fixed_Phi; /* azi (0N->90E) (rad) */ +extern float Telescope_x, Telescope_y; /* Telescope coordinates (cm) */ + extern int is_Random_Pointing; /* random pointing? */ extern float Random_Pointing_MaxDist; /* in metres */ @@ -49,8 +51,6 @@ 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) */ @@ -59,4 +59,5 @@ extern int ct_NPixels; /* number of pixels */ extern float ct_PixelWidth; /* pixel width (cm) */ + typedef struct Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c (revision 1431) @@ -10,4 +10,6 @@ extern char line[]; /* parsing buf. (init) */ +char axisdev_filename[256], reflectivity_filename[256]; + /* Prototypes */ extern void setall(long iseed1,long iseed2); /* rnds */ @@ -15,6 +17,6 @@ static void ReadCerfiles(FILE *parfile) -{ char *value_ptr = NULL; /* ptr at parm value */ - extern FILE *GetNextFile(char *cername); /* in main.c */ +{ char *value_ptr = NULL; /* ptr at parm value */ + extern FILE *GetNextFile(char *cername); /* in main.c */ filelist = parfile; @@ -22,11 +24,11 @@ if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL || (value_ptr=strtok(line, whites)) == NULL) - FatalError(FLST_NSPEC_FTL); + FatalError(FLST_NSPEC_FTL); else if (value_ptr[0] == '@') - { fclose(filelist); + { fclose(filelist); if ((filelist=fopen(value_ptr+1, "r")) == NULL) - FatalError(FLST_NFND__FTL, value_ptr+1); + FatalError(FLST_NFND__FTL, value_ptr+1); else if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) - FatalError(FLST_NSPEC_FTL); + FatalError(FLST_NSPEC_FTL); value_ptr = strtok(line, whites); } @@ -34,8 +36,8 @@ strcpy(cername, value_ptr); if ((value_ptr=strtok(NULL, whites)) == NULL) - { first_Event = 0; + { first_Event = 0; last_Event = 1000000; } else - { first_Event = atol(value_ptr); + { first_Event = atol(value_ptr); value_ptr = strtok(NULL, whites); last_Event = value_ptr ? atol(value_ptr) : 1000000; } @@ -44,18 +46,18 @@ if ((cerfile=fopen(cername, "r")) == NULL) - { Message(CERF_NFND__MSG, cername); - cerfile=GetNextFile(cername); } + { Message(CERF_NFND__MSG, cername); + cerfile=GetNextFile(cername); } /* If no valid cerfile is found then exit */ if (cerfile == NULL) - FatalError(CERF_NSPEC_FTL); + FatalError(CERF_NSPEC_FTL); /* Check boundaries */ if (first_Event > last_Event) - { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername); + { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername); first_Event = 0; last_Event = 1000000; } -} /* end of ReadCerfiles */ + } /* end of ReadCerfiles */ void ParmsSwitch(FILE *parfile) @@ -63,79 +65,93 @@ 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. */ + 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; + { 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 axisdev_file: + strcpy(axisdev_filename, value_ptr); + break; + case reflectivity_file: + strcpy(reflectivity_filename, value_ptr); + 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 cer_files: - ReadCerfiles(parfile); - switch_end = TRUE; - default: switch_end = TRUE; - break; }} + /* Added May 2002, AM: */ + case telescope_position: + Telescope_x = (float) atof(value_ptr); + value_ptr = strtok(NULL, whites); + if (value_ptr == NULL) + { Error(TEL_POS_ERR); + exit(-1);} + else + { Telescope_y = (float) atof(value_ptr); + Message(TEL_POS_MSG, Telescope_x, Telescope_y);} + 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 cer_files: + ReadCerfiles(parfile); + switch_end = TRUE; + default: switch_end = TRUE; + break; }} while (!switch_end); if (filelist == NULL) - FatalError(FLST_NSPEC_FTL); + FatalError(FLST_NSPEC_FTL); /* Set random seeds */ @@ -143,3 +159,3 @@ Message(SEED_SET___MSG, Seeds[0], Seeds[1]); -} /* end of ParmsSwitch */ + } /* end of ParmsSwitch */ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h (revision 1431) @@ -6,12 +6,13 @@ T(output_file), /* output file */ \ T(ct_file), /* file with the characteristics of the CT */ \ +T(axisdev_file), /* file with the single mirror spot deviations*/ \ +T(reflectivity_file), /* file with the mirror reflectivity */ \ 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(telescope_position),/* 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) */ @@ -39,4 +40,9 @@ #define FIXD_ENABL_MSG /* theta, phi */ \ "Using \"fixed_target\" mode with theta=%.2fdeg and phi=%.2fdeg.\n" +#define TEL_POS_ERR /* no parms */ \ + "Error while parsing \"telescope_position\" fields.\n" \ + " *** EXITING! \n" +#define TEL_POS_MSG /* theta, phi */ \ + "Using \"telescope position\" mode with x=%.1fcm and y=%.1fcm.\n" #define MAX__EVTS__MSG /* max. nr. of evts. */ \ "Processing at most %ld events.\n" @@ -59,3 +65,6 @@ " *** Please specify some *valid* filename after the cer_files directive.\n" + +extern char axisdev_filename[256], reflectivity_filename[256]; + #endif Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c (revision 1431) @@ -237,7 +237,14 @@ + 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]); + /* FIXED Lines below, May 2002, AM : formerly (up to V0.4) + * there was a confusion between telescope coordinates xCT and + * the original coordinates x. Thanks to T. Hengstebeck for + * reporting the bug. + */ + + c = 2*rCT[0]*rCT[2]*xCT[0]*xCT[2] + 2*rCT[1]*rCT[2]*xCT[1]*xCT[2] + - SQR(rCT[2])*SQR(xCT[0]) - SQR(rCT[2])*SQR(xCT[1]) + - SQR(rCT[0])*SQR(xCT[2]) - SQR(rCT[1])*SQR(xCT[2]); + if ( fabs(a) < 1.e-6 ) { @@ -583,5 +590,5 @@ Debug("@19 xcam-AD %f %f \n", xcam[0], xcam[1]); /* CBC */ - + /* ++ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c (revision 1421) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c (revision 1431) @@ -24,6 +24,9 @@ #include "header.h" +#define MAX(x,y) ((x)>(y)? (x) : (y)) +#define MIN(x,y) ((x)>(y)? (y) : (x)) + FILE *rflfile = NULL, /* reflector (output) file */ - *cerfile = NULL; /* current cerfile */ +*cerfile = NULL; /* current cerfile */ cphoton *CPhotons = NULL; /* cphoton array */ @@ -42,19 +45,26 @@ static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile); -FILE *chkf = NULL; /**********************/ -long myloop; /**********************/ - -void main(void) +FILE *chkf = NULL; +long myloop; + +Event_end* evt_end; + +main(void) { long event = 0L; /* event counter */ /* Read init & geometry parms, init files and vars */ init(NULL); - chkf = fopen("check", "w"); /**********************/ + + /* chkf = fopen("check", "w"); */ /* 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); } + { + if (ProcessEvent(cheadp, cerfile, rflfile)) + event++; + if (event % SHOW_ME == 0) + Log(INFO_EVENT_LOG, + cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal); + } /* Writing final flags */ @@ -64,5 +74,5 @@ Log(RFLF_CLOSE_LOG); fclose(rflfile); - fclose(chkf); /**********************/ + /* fclose(chkf); */ /* Clean memory and exit */ @@ -70,5 +80,5 @@ Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION)); -} /* end of main */ + } /* end of main */ static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile) @@ -82,23 +92,24 @@ extern void makeOmegaI(float theta, float phi); + char pp[256]; /* 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; + 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 */ + float wlen, /* photon wavelength */ + theta; /* photon zenith angle */ /* Reset counters and set fileptr */ @@ -111,24 +122,78 @@ /* Calculate OmegaCT matrices */ if (is_Fixed_Target) - { makeOmega(fixed_Theta, fixed_Phi); + { makeOmega(fixed_Theta, fixed_Phi); makeOmegaI(fixed_Theta, fixed_Phi); } else - { makeOmega(cheadp->Theta, cheadp->Phi); + { makeOmega(cheadp->Theta, cheadp->Phi); makeOmegaI(cheadp->Theta, cheadp->Phi); } memcpy( OmegaCT, Omega, 9*sizeof(float) ); memcpy( OmegaICT, OmegaI, 9*sizeof(float) ); + /* ADDED AM May 2002: now we take into account the telescope + * position chosen in the Corsika input card via the CERTEL + * option, which must be supplied also in the reflector input + * card with the option "telescope_position x y". Otherwise + * x = y = 0 is assumed, which is the standard mode of + * generation for MAGIC. + * Here we change rheadp so that the CorePos written to the + * .rfl file is referred to the telescope position. However, + * I believe that the original "CorePos" vector points + * from the core to the origin of coordinates of Corsika, and + * therefore after the subtraction of (Telescope_x,Telescope_y) + * the resulting vector CorePos points from the core to the + * telescope! + */ + + rheadp->CorePos[0][0] += Telescope_x; + rheadp->CorePos[1][0] += Telescope_y; + /* 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; phlongi_Nmax = evt_end->longi_Nmax; + rheadp->longi_t0 = evt_end->longi_t0; + rheadp->longi_tmax = evt_end->longi_tmax; + rheadp->longi_a = evt_end->longi_a; + rheadp->longi_b = evt_end->longi_b; + rheadp->longi_c = evt_end->longi_c; + rheadp->longi_chi2 = evt_end->longi_chi2; + break; + } + + for (ph=0; phCorePos[0][0]; - Photons[ph].y -= cheadp->CorePos[1][0]; + + /* TEMPORARY FIX, AM Nov 2001: we found that sometimes the value + stored in Photons[ph].w is not correct, and can result in a + wavelength beyond 600 nm, which makes the program crash later. + Now we force wlen to its expected range: + */ + + wlen = MIN(MAX(290.,wlen),600.); + + + /* ADDED AM May 2002: now we take into account the telescope + * position chosen in the Corsika input card via the CERTEL + * option, which must be supplied also in the reflector input + * card with the option "telescope_position x y". Otherwise + * x = y = 0 is assumed, which is the standard mode of + * generation for MAGIC. + */ + + Photons[ph].x -= cheadp->CorePos[0][0]+Telescope_x; + Photons[ph].y -= cheadp->CorePos[1][0]+Telescope_y; /* Increment number of photons */ @@ -137,22 +202,24 @@ /* Calculate some quantities */ theta = (float) acos(sqrt( - 1.f - Photons[ph].u*Photons[ph].u - - Photons[ph].v*Photons[ph].v)); + MAX(0., 1.f - Photons[ph].u*Photons[ph].u + - Photons[ph].v*Photons[ph].v))); /* Check absorption */ + + if (absorption(wlen, Photons[ph].h, theta)) - absphs++; + absphs++; /* Check reflection */ else if (0 != (ref_type = ph2cph(&Photons[ph], &CPhotons[cphs]))) - refphs[ref_type-1]++; - else /* Photon passed */ - { - Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph, - Photons[ph].x,Photons[ph].y, - Photons[ph].u,Photons[ph].v,theta); - Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x, - CPhotons[cphs].y); + refphs[ref_type-1]++; + else /* Photon passed */ + { + Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph, + Photons[ph].x,Photons[ph].y, + Photons[ph].u,Photons[ph].v,theta); + Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x, + CPhotons[cphs].y); @@ -162,9 +229,9 @@ /* Update cphs */ - if (++cphs == NR_OF_CPHOTONS) /* Overflow */ - { + if (++cphs == NR_OF_CPHOTONS) /* Overflow */ + { /* if it is the first o/f open a tempfile */ if (overflow++ == 0) - { Log("Spooling... "); + { Log("Spooling... "); tmpf = tmpfile(); if (tmpf == NULL) FatalError(TEMP_ERROR_FTL); } @@ -176,15 +243,15 @@ cphs = 0; - } /* if overflow */ - } /* else (=Photon passed) */ - } /* end of loop inside datablock */ - } /* end of loop on datablocks */ + } /* 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 */ + else /* no error: write the new event */ + + { /* Write "start of event" flag */ fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile); @@ -202,18 +269,20 @@ rheadp->CPhotons = (long) overflow * NR_OF_CPHOTONS + cphs; fwrite(rheadp, sizeof(RflHeader), 1, rflfile); - for (myloop=0; myloop