Changeset 1431 for trunk/MagicSoft
- Timestamp:
- 07/24/02 15:36:10 (22 years ago)
- Location:
- trunk/MagicSoft/Simulation/Detector/ReflectorII
- Files:
-
- 13 edited
-
Makefile (modified) (9 diffs)
-
attenu.f (modified) (3 diffs)
-
geometry.c (modified) (7 diffs)
-
geometry.h (modified) (3 diffs)
-
header.c (modified) (1 diff)
-
header.h (modified) (5 diffs)
-
init.c (modified) (1 diff)
-
init.h (modified) (4 diffs)
-
parms.c (modified) (7 diffs)
-
parms.h (modified) (3 diffs)
-
ph2cph.c (modified) (2 diffs)
-
reflector.c (modified) (15 diffs)
-
version.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile
r774 r1431 19 19 # 20 20 # $RCSfile: Makefile,v $ 21 # $Revision: 1. 2$22 # $Author: b lanch$23 # $Date: 200 1-04-26 07:57:41$21 # $Revision: 1.3 $ 22 # $Author: bigongia $ 23 # $Date: 2002-07-24 14:35:45 $ 24 24 # 25 25 ################################################################## … … 53 53 # compilation and linking flags 54 54 55 CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} ${DEBUG} 55 #CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} ${DEBUG} 56 CXXFLAGS = -D__${SYSTEM}__ ${INCLUDES} ${OPTIM} 56 57 CFLAGS = ${CXXFLAGS} 57 #CFLAGS = -std1 -warnprotos -msg_enable level6 -msg_disable returnchecks58 58 FFLAGS = ${CXXFLAGS} 59 LIBS = ${RANLIB} ${FORLIBS}60 59 LIBS = ${CERNLIB} ${RANLIB} ${FORLIBS} 61 60 … … 180 179 diag.o: /usr/include/bits/types.h /usr/include/libio.h 181 180 diag.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h 182 diag.o: /usr/include/stdlib.h /usr/include/sys/types.h /usr/include/time.h 183 diag.o: /usr/include/endian.h /usr/include/bits/endian.h 184 diag.o: /usr/include/sys/select.h /usr/include/bits/select.h 185 diag.o: /usr/include/bits/sigset.h /usr/include/sys/sysmacros.h 186 diag.o: /usr/include/alloca.h version.h diag.h 181 diag.o: /usr/include/stdlib.h version.h diag.h 187 182 init.o: /usr/include/stdio.h /usr/include/features.h /usr/include/sys/cdefs.h 188 183 init.o: /usr/include/gnu/stubs.h … … 191 186 init.o: /usr/include/bits/types.h /usr/include/libio.h 192 187 init.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h 193 init.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/sys/types.h 194 init.o: /usr/include/time.h /usr/include/endian.h /usr/include/bits/endian.h 195 init.o: /usr/include/sys/select.h /usr/include/bits/select.h 196 init.o: /usr/include/bits/sigset.h /usr/include/sys/sysmacros.h 197 init.o: /usr/include/alloca.h /usr/include/math.h 188 init.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/math.h 198 189 init.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h 199 init.o: /usr/include/bits/mathcalls.h 200 init.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h 201 init.o: version.h diag.h init.h 190 init.o: /usr/include/bits/mathcalls.h version.h diag.h init.h 202 191 parms.o: /usr/include/stdio.h /usr/include/features.h 203 192 parms.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h … … 206 195 parms.o: /usr/include/bits/types.h /usr/include/libio.h 207 196 parms.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h 208 parms.o: /usr/include/string.h /usr/include/ math.h197 parms.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/math.h 209 198 parms.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h 210 parms.o: /usr/include/bits/mathcalls.h 211 parms.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h 212 parms.o: diag.h parms.h init.h 199 parms.o: /usr/include/bits/mathcalls.h diag.h parms.h init.h 213 200 geometry.o: /usr/include/stdio.h /usr/include/features.h 214 201 geometry.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h … … 217 204 geometry.o: /usr/include/bits/types.h /usr/include/libio.h 218 205 geometry.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h 219 geometry.o: /usr/include/string.h /usr/include/ math.h206 geometry.o: /usr/include/string.h /usr/include/stdlib.h /usr/include/math.h 220 207 geometry.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h 221 geometry.o: /usr/include/bits/mathcalls.h 222 geometry.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h 223 geometry.o: diag.h geometry.h init.h 208 geometry.o: /usr/include/bits/mathcalls.h diag.h geometry.h init.h 224 209 atm.o: /usr/include/stdio.h /usr/include/features.h /usr/include/sys/cdefs.h 225 210 atm.o: /usr/include/gnu/stubs.h … … 229 214 atm.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h 230 215 atm.o: /usr/include/string.h /usr/include/math.h /usr/include/bits/huge_val.h 231 atm.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h 232 atm.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h diag.h 233 atm.o: atm.h init.h 216 atm.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h diag.h atm.h 217 atm.o: init.h 234 218 ph2cph.o: /usr/include/stdio.h /usr/include/features.h 235 219 ph2cph.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h … … 239 223 ph2cph.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h 240 224 ph2cph.o: /usr/include/math.h /usr/include/bits/huge_val.h 241 ph2cph.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h 242 ph2cph.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h 243 ph2cph.o: diag.h init.h lagrange.h 225 ph2cph.o: /usr/include/bits/mathdef.h /usr/include/bits/mathcalls.h diag.h 226 ph2cph.o: init.h lagrange.h 244 227 header.o: /usr/include/string.h /usr/include/features.h 245 228 header.o: /usr/include/sys/cdefs.h /usr/include/gnu/stubs.h … … 252 235 reflector.o: /usr/include/bits/types.h /usr/include/libio.h 253 236 reflector.o: /usr/include/_G_config.h /usr/include/bits/stdio_lim.h 254 reflector.o: /usr/include/stdlib.h /usr/include/sys/types.h 255 reflector.o: /usr/include/time.h /usr/include/endian.h 256 reflector.o: /usr/include/bits/endian.h /usr/include/sys/select.h 257 reflector.o: /usr/include/bits/select.h /usr/include/bits/sigset.h 258 reflector.o: /usr/include/sys/sysmacros.h /usr/include/alloca.h 259 reflector.o: /usr/include/string.h /usr/include/math.h 237 reflector.o: /usr/include/stdlib.h /usr/include/string.h /usr/include/math.h 260 238 reflector.o: /usr/include/bits/huge_val.h /usr/include/bits/mathdef.h 261 reflector.o: /usr/include/bits/mathcalls.h 262 reflector.o: /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include/float.h 263 reflector.o: version.h diag.h init.h header.h 239 reflector.o: /usr/include/bits/mathcalls.h version.h diag.h init.h header.h -
trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f
r725 r1431 360 360 x4 = sqrt((2. * h + rt) / (2. * hscale)) 361 361 362 c-- AM Dec 2001, to avoid crash! A few photons seem to be "corrupted" 363 c-- (have absurd value) in a cer file... 364 365 if (abs(x3-x4) .lt. 1.e-10) then 366 tr_atmos = -1. 367 RETURN 368 endif 369 370 362 371 e1 = derfc(x1) 363 372 e2 = derfc(x2) … … 366 375 367 376 m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4)) 377 368 378 369 379 ********************************************************************** … … 442 452 443 453 CON_OZ = 2 454 444 455 2001 IF (LONG(CON_OZ).GT.wavelength) THEN 445 A = (OZ_ABI(ROW,CON_OZ)-OZ_ABI(ROW,CON_OZ-1))/456 A = (OZ_ABI(ROW,CON_OZ)-OZ_ABI(ROW,CON_OZ-1))/ 446 457 + (LONG(CON_OZ) - LONG(CON_OZ-1)) 447 458 B = OZ_ABI(ROW,CON_OZ) - (A*LONG(CON_OZ)) -
trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.c
r923 r1431 8 8 9 9 extern char line[]; /* parsing buf. (init) */ 10 extern char axisdev_filename[256], reflectivity_filename[256]; 11 12 float mean_refl; /* Mirror mean reflectivity 270-610 nm. 13 * AM June 2002. 14 */ 10 15 11 16 float ct_Focal_mean; /* focal dist. (mean) (cm) */ 12 float ct_Focal_std; /* focal dist. (std) (cm) */13 17 float ct_PSpread_mean; /* pt. spread fn. (mean) (cm) */ 14 float ct_PSpread_std; /* pt. spread fn. (std) (cm) */15 float ct_Adjustment_std; /* adjustment dev. (std) (cm) */16 18 float ct_BlackSpot_rad; /* black spot radius (cm) */ 17 19 float ct_RMirror; /* rad. of single mirror (cm) */ … … 43 45 Log(MIRR_ALLOC_LOG, ct_NMirrors); 44 46 Log(MIRR_TABLE_LOG); 45 if (ct_BinaryData) 46 { Log(BINF_OPEN__LOG, ct_BinaryName); 47 fread(ct_data, sizeof(mirror), ct_NMirrors, ct_BinaryData); 48 fclose(ct_BinaryData); } 49 else 50 { /* read ASCII data */ 51 Log(READ_ASCII_LOG); 52 for (i=0; i<ct_NMirrors; i++) 53 { if (12 != fscanf(geofile, "%d %f %f %f %f %f %f %f %f %f %f %f", 54 &ct_data[i].i, &ct_data[i].f, 55 &ct_data[i].sx, &ct_data[i].sy, 56 &ct_data[i].x, &ct_data[i].y, &ct_data[i].z, 57 &ct_data[i].theta, &ct_data[i].phi, 58 &ct_data[i].xn, &ct_data[i].yn, &ct_data[i].zn)) 59 break; 60 Log("[%d]", ct_data[i].i); } 61 Log("%c", '\n'); 62 if (i < ct_NMirrors) 63 FatalError(MIRR_FEW___FTL, i); 64 /* Data Ok: save binary data for next time */ 65 if ((ct_BinaryData=fopen(ct_BinaryName, "w"))==NULL) 66 Log(BINF_ERROR_LOG, ct_BinaryName); 67 else 68 { Log(BINF_WRITE_LOG, ct_BinaryName); 69 fwrite(ct_data, sizeof(mirror), ct_NMirrors, ct_BinaryData); 70 fclose(ct_BinaryData); } 71 } /* end of if: reading ASCII data */ 47 48 Log(READ_ASCII_LOG); 49 for (i=0; i<ct_NMirrors; i++) 50 { 51 if (12 != fscanf(geofile, "%d %f %f %f %f %f %f %f %f %f %f %f", 52 &ct_data[i].i, &ct_data[i].f, 53 &ct_data[i].sx, &ct_data[i].sy, 54 &ct_data[i].x, &ct_data[i].y, &ct_data[i].z, 55 &ct_data[i].theta, &ct_data[i].phi, 56 &ct_data[i].xn, &ct_data[i].yn, &ct_data[i].zn)) 57 break; 58 } 59 if (i < ct_NMirrors) 60 FatalError(MIRR_FEW___FTL, i); 61 72 62 } /* end of ReadMirrorTable */ 73 63 74 64 static void ReadReflectivity(char *datname) 75 { FILE *datfile = fopen(datname, "r"); 65 { 66 FILE *datfile = fopen(datname, "r"); 76 67 int current = 0; 68 69 mean_refl = 0.; 77 70 78 71 if (datfile == NULL) 79 72 FatalError(RFLF_ERROR_FTL, datname); 73 else 74 printf("Reading file %s\n", datname); 75 80 76 while (fgets(line, LINE_MAX_LENGTH, datfile)) 81 { if (line[0] == '#') continue; 82 if (nReflectivity == 0) 83 { nReflectivity = atoi(line); 84 if (nReflectivity) 85 { if ((Reflectivity[0] = 86 (float *) malloc(sizeof(float) * nReflectivity)) == NULL 87 || (Reflectivity[1] = 88 (float *) malloc(sizeof(float) * nReflectivity)) == NULL) 89 FatalError(REFL_ALLOC_FTL, nReflectivity); }} 90 else if (2 == sscanf(line, "%f %f", 91 &Reflectivity[0][current], &Reflectivity[1][current])); 92 { current++; 93 if (current >= nReflectivity) break; }} 77 { 78 if (line[0] == '#') continue; 79 80 if (nReflectivity == 0) 81 { 82 nReflectivity = atoi(line); 83 if (nReflectivity) 84 { 85 if ((Reflectivity[0] = 86 (float *) malloc(sizeof(float) * nReflectivity)) == NULL || 87 (Reflectivity[1] = 88 (float *) malloc(sizeof(float) * nReflectivity)) == NULL) 89 FatalError(REFL_ALLOC_FTL, nReflectivity); 90 } 91 } 92 else if (2 == sscanf(line, "%f %f", &Reflectivity[0][current], 93 &Reflectivity[1][current])) 94 { 95 // Added June 2002, AM: 96 mean_refl += Reflectivity[1][current]; 97 current++; 98 if (current >= nReflectivity) break; 99 } 100 } 94 101 fclose(datfile); 95 102 96 103 nReflectivity = current; 104 if (current > 0) 105 mean_refl /= (float) current; 106 97 107 } /* end of ReadReflectivity */ 98 108 … … 103 113 if (datfile == NULL) 104 114 FatalError(AXIS_ERROR_FTL, datname); 115 else 116 printf("Reading file %s\n", axisdev_filename); 117 105 118 if ((AxisDeviation[0]= 106 119 (float *) malloc(sizeof(float) * ct_NMirrors)) == NULL … … 149 162 break; 150 163 case focal_distance: 151 Log(LOG__FLOAT_LOG, "focal distance (average )",164 Log(LOG__FLOAT_LOG, "focal distance (average, cm)", 152 165 ct_Focal_mean = (float) atof(value_ptr)); 153 166 break; 154 case focal_std: 155 Log(LOG__FLOAT_LOG, "focal distance (std. dev.)", 156 ct_Focal_std = (float) atof(value_ptr)); 167 case focal_std: /* not implemented. */ 157 168 break; 158 169 case point_spread: 159 Log(LOG__FLOAT_LOG, "point spread fn. (average)",170 Log(LOG__FLOAT_LOG, "point spread fn. sigma (average, cm)", 160 171 ct_PSpread_mean = (float) atof(value_ptr)); 161 172 break; 162 case point_std: 163 Log(LOG__FLOAT_LOG, "point spread fn. (std. dev.)", 164 ct_PSpread_std = (float) atof(value_ptr)); 165 break; 166 case adjustment_dev: 167 Log(LOG__FLOAT_LOG, "adjustment dev. (std. dev.)", 168 ct_Adjustment_std = (float) atof(value_ptr)); 173 case point_std: /* not implemented */ 174 break; 175 case adjustment_dev: /* not implemented */ 169 176 break; 170 177 case black_spot: … … 172 179 ct_BlackSpot_rad = (float) atof(value_ptr)); 173 180 break; 181 case n_mirrors: 182 Log(LOG__INT___LOG, "number of mirrors", 183 ct_NMirrors = atoi(value_ptr)); 184 break; 174 185 case r_mirror: 175 186 Log(LOG__FLOAT_LOG, "single mirror radius (cm)", 176 187 ct_RMirror = (float) atof(value_ptr)); 177 188 break; 178 case n_mirrors:179 Log(LOG__INT___LOG, "number of mirrors",180 ct_NMirrors = atoi(value_ptr));181 break;182 189 case camera_width: 183 Log(LOG__FLOAT_LOG, "camera width(cm)",190 Log(LOG__FLOAT_LOG, "camera radius (cm)", 184 191 ct_CameraWidth = (float) atof(value_ptr)); 185 192 break; … … 191 198 Log(LOG__FLOAT_LOG, "pixel width (cm)", 192 199 ct_PixelWidth = (float) atof(value_ptr)); 200 break; 201 case n_centralpixels: 202 /* this parameter is for camera, not for reflector */ 203 break; 204 case n_gappixels: 205 /* this parameter is for camera, not for reflector */ 193 206 break; 194 207 case refl_file: … … 208 221 fclose(geofile); 209 222 210 if (Reflectivity[0] == NULL) ReadReflectivity(REFLECTIVITY_FILE); 211 if (AxisDeviation[0]== NULL) ReadAxisDev(AXISDEVIATION_FILE); 223 if (strlen(reflectivity_filename) == 0) 224 strcpy(reflectivity_filename, REFLECTIVITY_FILE); 225 if (Reflectivity[0] == NULL) ReadReflectivity(reflectivity_filename); 226 227 if (strlen(axisdev_filename) == 0) 228 strcpy(axisdev_filename, AXISDEVIATION_FILE); 229 if (AxisDeviation[0]== NULL) ReadAxisDev(axisdev_filename); 230 231 212 232 ReadFocals(); 213 233 } /* end of GeometrySwitch */ 234 -
trunk/MagicSoft/Simulation/Detector/ReflectorII/geometry.h
r725 r1431 17 17 T(refl_file), /* path of file containing refl. data */ \ 18 18 T(axisdev_file), /* path of file containing axis dev. data */ \ 19 T(define_mirrors) /* this entry is followed by the def. of pixels */ 20 19 T(define_mirrors), /* this entry is followed by the def. of pixels */ \ 20 T(n_centralpixels), /* this token is not for reflector but for camera */ \ 21 T(n_gappixels) /* same comment as for previous token. */ 22 21 23 #define T(x) x /* define T() as the name as it is */ 22 24 enum { ITEM_LIST }; … … 48 50 "Located and opened file \"%s\" containing mirror data.\n" 49 51 #define READ_ASCII_LOG /* no parms */ \ 50 " Start reading ASCII data for mirror nr:\n"52 "Reading ASCII data for mirrors.\n" 51 53 #define MIRR_FEW___FTL /* mirrors read */ \ 52 54 "Not enough mirror data: only %d mirrors read.\n" 53 #define BINF_ERROR_LOG /* bin. data filename */ \54 "Cannot write mirror data on binary file:\n %s\n" \55 " *** Skipping operation and proceeding.\n"56 #define BINF_WRITE_LOG /* bin. data filename */ \57 "Binary data written in file \"%s\".\n"58 55 #define RFLF_ERROR_FTL /* reflectivity fname */ \ 59 56 "Cannot find file \"%s\" containing data on reflectivity.\n" … … 82 79 */ 83 80 81 extern float mean_refl; /* mean mirror reflectivity AM June 2002 */ 82 84 83 #define RANDOM_POINTING_MAX_SEPARATION 0.104719755119660 85 84 -
trunk/MagicSoft/Simulation/Detector/ReflectorII/header.c
r725 r1431 9 9 static CerHeader chead; CerHeader *cheadp = &chead; 10 10 11 extern float fixed_Phi, fixed_Theta; 12 extern int ct_NMirrors; 13 extern float mean_refl; 14 11 15 void TranslateHeader(RflHeader *r, CerHeader *c) 12 16 { 13 r->EvtNumber = c->EvtNumber; 14 r->PrimaryID = c->PrimaryID; 15 r->Etotal = c->Etotal; 16 r->Thick0 = c->Thick0; 17 r->FirstTarget = c->FirstTarget; 18 r->zFirstInt = c->zFirstInt; 19 r->Theta = c->Theta; 20 r->Phi = c->Phi; 21 r->NumRndSeq = c->NumRndSeq; 22 r->RunNumber = c->RunNumber; 23 r->DateRun = c->DateRun; 24 r->VersionPGM = c->VersionPGM; 25 r->NumObsLev = c->NumObsLev; 26 r->SlopeSpec = c->SlopeSpec; 27 r->ELowLim = c->ELowLim; 28 r->EUppLim = c->EUppLim; 29 r->ThetaMin = c->ThetaMin; 30 r->ThetaMax = c->ThetaMax; 31 r->PhiMin = c->PhiMin; 32 r->PhiMax = c->PhiMax; 33 r->CWaveLower = c->CWaveLower; 34 r->CWaveUpper = c->CWaveUpper; 17 r->EvtNumber = c->EvtNumber; 18 r->PrimaryID = c->PrimaryID; 19 r->Etotal = c->Etotal; 20 r->Thick0 = c->Thick0; 21 r->FirstTarget = c->FirstTarget; 22 r->zFirstInt = c->zFirstInt; 23 r->Theta = c->Theta; 24 r->Phi = c->Phi; 25 r->NumRndSeq = c->NumRndSeq; 26 r->RunNumber = c->RunNumber; 27 r->DateRun = c->DateRun; 28 r->Corsika_version= c->Corsika_version; 29 r->NumObsLev = c->NumObsLev; 30 r->HeightLev = c->HeightLev[0]; 35 31 36 memcpy(r->p, c->p, 3*sizeof(float));37 memcpy(r->RndData, c->RndData, 30*sizeof(float));38 memcpy(r->HeightLev, c->HeightLev, 10*sizeof(float));39 memcpy(r->CorePos, c->CorePos, 40*sizeof(float));40 32 41 r->deviationPhi = r->deviationTheta = r->Trigger = 0.f; 33 r->SlopeSpec = c->SlopeSpec; 34 r->ELowLim = c->ELowLim; 35 r->EUppLim = c->EUppLim; 36 r->ThetaMin = c->ThetaMin; 37 r->ThetaMax = c->ThetaMax; 38 r->PhiMin = c->PhiMin; 39 r->PhiMax = c->PhiMax; 40 r->CWaveLower = c->CWaveLower; 41 r->CWaveUpper = c->CWaveUpper; 42 42 43 } /* end of TranslateHeader */ 43 memcpy(r->p, c->p, 3*sizeof(float)); 44 memcpy(r->RndData, c->RndData, 30*sizeof(float)); 45 memcpy(r->CorePos, c->CorePos, 40*sizeof(float)); 46 47 /* Next 4 variables added in June 2002, AM */ 48 r->telescopePhi = fixed_Phi; 49 r->telescopeTheta = fixed_Theta; 50 r->num_mirrors = ct_NMirrors; 51 r->mean_reflectivity = mean_refl; 52 53 if (c->Corsika_version >= 6.) /* Viewcone option implemented only in c6xx */ 54 { 55 if (c->viewcone_angles[0] > 0.) 56 { 57 printf ("ERROR: Input cer file was run using Corsika option\n"); 58 printf (" VIEWCONE VUECON(1) VUECON(2) with VUECON(1) > 0\n"); 59 printf (" The reflector program only supports VUECON(1) = 0\n\n"); 60 exit(-1); 61 } 62 else 63 r->ViewConeRadius = c->viewcone_angles[1]; /* degrees */ 64 } 65 else 66 r->ViewConeRadius = 0.; 67 68 } /* end of TranslateHeader */ 69 70 -
trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h
r725 r1431 19 19 float RunNumber; 20 20 float DateRun; 21 float VersionPGM;21 float Corsika_version; 22 22 23 23 float NumObsLev; … … 72 72 float CorePos[2][20]; 73 73 74 float dmmy1; 75 float SpinTheta; 76 float SpinPhi; 77 float dmmy2[132]; 74 float dmmy1[14]; 75 float viewcone_angles[2]; /* (degrees) Inner and outer angles in 76 * Corsika's VIEWCONE option 77 */ 78 float dmmy2[119]; 78 79 } CerHeader; 79 80 … … 94 95 float RunNumber; 95 96 float DateRun; 96 float VersionPGM;97 float Corsika_version; 97 98 98 float NumObsLev; 99 float HeightLev[10]; 99 float NumObsLev; /* Should be 1 for MAGIC simulation */ 100 float HeightLev; /* Observation Level */ 101 102 /* Changed meaning of next 9 variables. June 2002, A.Moralejo: */ 103 float num_mirrors; 104 float mean_reflectivity; 105 float longi_Nmax; 106 float longi_t0; 107 float longi_tmax; 108 float longi_a; 109 float longi_b; 110 float longi_c; 111 float longi_chi2; 100 112 101 113 float SlopeSpec; … … 114 126 float TimeLast; 115 127 116 float deviationPhi; 117 float deviationTheta; 128 /* AM, 23/05/2002: Changed meaning of following 129 * three variables (which were unused before): 130 */ 131 float telescopePhi; /* rad */ 132 float telescopeTheta; /* rad */ 118 133 119 float Trigger; 134 float ViewConeRadius; /* Degrees. 135 * Radius of "view cone" when the primaries' 136 * directions generated by Corsika lie within 137 * a cone around a fixed direction. This is 138 * only possible with Corsika>6 versions. In 139 * that case, PhiMin=PhiMax and 140 * ThetaMin=ThetaMax (also in this header) 141 * indicate the axis of this cone. 142 * If ViewConeRadius==0, it means that 143 * the VIEWCONE option was not used. 144 */ 145 120 146 121 147 float CORSIKAPhs; /* Original photons written by Corsika */ … … 124 150 float OutOfMirrPhs; /* Photons outside the mirror */ 125 151 float BlackSpotPhs; /* Photons lost in the "black spot" */ 126 float OutOfChamPhs; /* Photons outside the c hamber*/127 float CPhotons; /* Photons reaching the c hamber*/152 float OutOfChamPhs; /* Photons outside the camera */ 153 float CPhotons; /* Photons reaching the camera */ 128 154 } RflHeader; 155 156 157 typedef struct /* EVTE from cerfile */ 158 { 159 float dummy1[255]; 160 float longi_Nmax; 161 float longi_t0; 162 float longi_tmax; 163 float longi_a; 164 float longi_b; 165 float longi_c; 166 float longi_chi2; 167 float dummy2[11]; 168 } Event_end; 129 169 130 170 void TranslateHeader(RflHeader *r, CerHeader *c); -
trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c
r725 r1431 20 20 float fixed_Theta, /* zenith angle (rad) */ 21 21 fixed_Phi; /* azi (0N->90E) (rad) */ 22 23 float Telescope_x = 0.; 24 float Telescope_y = 0.; /* Telescope coordinates (cm) */ 22 25 23 26 int is_Random_Pointing=FALSE; /* random pointing? */ -
trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h
r774 r1431 4 4 /* Constant definitions */ 5 5 #define LINE_MAX_LENGTH 128 /* parsing buffer sz. */ 6 #define NR_OF_CPHOTONS 65535/* CPhotons sz. */6 #define NR_OF_CPHOTONS 65535 /* CPhotons sz. */ 7 7 #define PH_IN_DATABLOCK 39 /* Photons in datablk */ 8 8 #define SHOW_ME 500 /* How many evts among two logs */ … … 41 41 extern float fixed_Phi; /* azi (0N->90E) (rad) */ 42 42 43 extern float Telescope_x, Telescope_y; /* Telescope coordinates (cm) */ 44 43 45 extern int is_Random_Pointing; /* random pointing? */ 44 46 extern float Random_Pointing_MaxDist; /* in metres */ … … 49 51 extern char ct_BinaryName[]; /* binary data filename */ 50 52 extern float ct_Focal_mean; /* focal dist. (mean) (cm) */ 51 extern float ct_Focal_std; /* focal dist. (std) (cm) */52 53 extern float ct_PSpread_mean; /* pt. spread fn. (mean) (cm) */ 53 54 extern float ct_PSpread_std; /* pt. spread fn. (std) (cm) */ 54 extern float ct_Adjustment_std; /* adjustment dev. (std) (cm) */55 55 extern float ct_BlackSpot_rad; /* black spot radius (cm) */ 56 56 extern float ct_RMirror; /* rad. of single mirror (cm) */ … … 59 59 extern int ct_NPixels; /* number of pixels */ 60 60 extern float ct_PixelWidth; /* pixel width (cm) */ 61 61 62 62 63 typedef struct -
trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c
r923 r1431 10 10 extern char line[]; /* parsing buf. (init) */ 11 11 12 char axisdev_filename[256], reflectivity_filename[256]; 13 12 14 /* Prototypes */ 13 15 extern void setall(long iseed1,long iseed2); /* rnds */ … … 15 17 16 18 static void ReadCerfiles(FILE *parfile) 17 { char *value_ptr = NULL; 18 extern FILE *GetNextFile(char *cername); 19 { char *value_ptr = NULL; /* ptr at parm value */ 20 extern FILE *GetNextFile(char *cername); /* in main.c */ 19 21 20 22 filelist = parfile; … … 22 24 if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL || 23 25 (value_ptr=strtok(line, whites)) == NULL) 24 26 FatalError(FLST_NSPEC_FTL); 25 27 else if (value_ptr[0] == '@') 26 { fclose(filelist);28 { fclose(filelist); 27 29 if ((filelist=fopen(value_ptr+1, "r")) == NULL) 28 30 FatalError(FLST_NFND__FTL, value_ptr+1); 29 31 else if (fgets(line, LINE_MAX_LENGTH, filelist) == NULL) 30 32 FatalError(FLST_NSPEC_FTL); 31 33 value_ptr = strtok(line, whites); } 32 34 … … 34 36 strcpy(cername, value_ptr); 35 37 if ((value_ptr=strtok(NULL, whites)) == NULL) 36 { first_Event = 0;38 { first_Event = 0; 37 39 last_Event = 1000000; } 38 40 else 39 { first_Event = atol(value_ptr);41 { first_Event = atol(value_ptr); 40 42 value_ptr = strtok(NULL, whites); 41 43 last_Event = value_ptr ? atol(value_ptr) : 1000000; } … … 44 46 45 47 if ((cerfile=fopen(cername, "r")) == NULL) 46 { Message(CERF_NFND__MSG, cername);47 cerfile=GetNextFile(cername); }48 { Message(CERF_NFND__MSG, cername); 49 cerfile=GetNextFile(cername); } 48 50 49 51 /* If no valid cerfile is found then exit */ 50 52 if (cerfile == NULL) 51 53 FatalError(CERF_NSPEC_FTL); 52 54 53 55 /* Check boundaries */ 54 56 if (first_Event > last_Event) 55 { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername);57 { Error(EVTN_WRONG_ERR, first_Event, last_Event, cername); 56 58 first_Event = 0; 57 59 last_Event = 1000000; } 58 60 59 }/* end of ReadCerfiles */61 } /* end of ReadCerfiles */ 60 62 61 63 void ParmsSwitch(FILE *parfile) … … 63 65 int switch_end = FALSE; /* bool to exit loop */ 64 66 extern FILE *geofile; /* geo file (init) */ 65 extern void SetVerbose(int vlevel); 66 extern void SetAtmModel(char *model); 67 extern int ParseLine(FILE *parfile, 68 const char *token_list[], 69 int tokens, 70 char **value_ptr); 67 extern void SetVerbose(int vlevel); /* from diag.c */ 68 extern void SetAtmModel(char *model); /* from atm.c */ 69 extern int ParseLine(FILE *parfile, /* FILE with parms */ 70 const char *token_list[], /* array w/tokens */ 71 int tokens, /* nr of tokens */ 72 char **value_ptr); /* ptr->parm val. */ 71 73 72 74 do 73 { switch(ParseLine(parfile, parms, ARRAY_SZ(parms), &value_ptr)) 74 { case output_file: 75 if ((rflfile=fopen(value_ptr, "w+")) == NULL) 76 FatalError(OUTF_ERROR_FTL, value_ptr); 77 Message(OUTF_OPEN__MSG, value_ptr); 78 break; 79 case ct_file: 80 if ((geofile=fopen(value_ptr, "r")) == NULL) 81 FatalError(GEOF_ERROR_FTL, value_ptr); 82 Message(GEOF_OPEN__MSG, value_ptr); 83 strcat(strcpy(ct_BinaryName, value_ptr), ".mirr"); 84 ct_BinaryData = fopen(ct_BinaryName, "r"); 85 break; 86 case atm_model: 87 SetAtmModel(value_ptr); 88 break; 89 case verbose_level: 90 SetVerbose(atoi(value_ptr)); 91 break; 92 case fixed_target: 93 is_Fixed_Target = TRUE; 94 fixed_Theta = (float) atof(value_ptr); 95 value_ptr = strtok(NULL, whites); 96 if (value_ptr == NULL) 97 { Error(FIXD_TARGT_ERR); 98 is_Fixed_Target = FALSE; } 99 else 100 { fixed_Phi = (float) atof(value_ptr); 101 Message(FIXD_ENABL_MSG, fixed_Theta, fixed_Phi); 102 fixed_Theta *= (float) (M_PI/180.); 103 fixed_Phi *= (float) (M_PI/180.); } 104 break; 105 case max_events: 106 Message(MAX__EVTS__MSG, max_Events=atol(value_ptr)); 107 break; 108 case energy_cuts: 109 low_Ecut = (float) atof(value_ptr); 110 value_ptr = strtok(NULL, whites); 111 if (value_ptr == NULL) 112 { Error(ENRG_LIMIT_ERR); 113 low_Ecut = 0.; } 114 else 115 { high_Ecut = (float) atof(value_ptr); 116 Message(ENRG_CUTS__MSG, low_Ecut, high_Ecut); } 117 break; 118 case seeds: 119 Seeds[0] = atol(value_ptr); 120 value_ptr = strtok(NULL, whites); 121 if (value_ptr) Seeds[1] = atol(value_ptr); 122 else 123 { Error(SEED_ERROR_ERR); 124 Seeds[0] = 3141592L; } 125 break; 126 case random_pointing: 127 case repeat_random: 128 /********************************/ 129 break; 75 { switch(ParseLine(parfile, parms, ARRAY_SZ(parms), &value_ptr)) 76 { case output_file: 77 if ((rflfile=fopen(value_ptr, "w+")) == NULL) 78 FatalError(OUTF_ERROR_FTL, value_ptr); 79 Message(OUTF_OPEN__MSG, value_ptr); 80 break; 81 case ct_file: 82 if ((geofile=fopen(value_ptr, "r")) == NULL) 83 FatalError(GEOF_ERROR_FTL, value_ptr); 84 Message(GEOF_OPEN__MSG, value_ptr); 85 strcat(strcpy(ct_BinaryName, value_ptr), ".mirr"); 86 ct_BinaryData = fopen(ct_BinaryName, "r"); 87 break; 88 case axisdev_file: 89 strcpy(axisdev_filename, value_ptr); 90 break; 91 case reflectivity_file: 92 strcpy(reflectivity_filename, value_ptr); 93 break; 94 case atm_model: 95 SetAtmModel(value_ptr); 96 break; 97 case verbose_level: 98 SetVerbose(atoi(value_ptr)); 99 break; 100 case fixed_target: 101 is_Fixed_Target = TRUE; 102 fixed_Theta = (float) atof(value_ptr); 103 value_ptr = strtok(NULL, whites); 104 if (value_ptr == NULL) 105 { Error(FIXD_TARGT_ERR); 106 is_Fixed_Target = FALSE; } 107 else 108 { fixed_Phi = (float) atof(value_ptr); 109 Message(FIXD_ENABL_MSG, fixed_Theta, fixed_Phi); 110 fixed_Theta *= (float) (M_PI/180.); 111 fixed_Phi *= (float) (M_PI/180.); } 112 break; 130 113 131 case cer_files: 132 ReadCerfiles(parfile); 133 switch_end = TRUE; 134 default: switch_end = TRUE; 135 break; }} 114 /* Added May 2002, AM: */ 115 case telescope_position: 116 Telescope_x = (float) atof(value_ptr); 117 value_ptr = strtok(NULL, whites); 118 if (value_ptr == NULL) 119 { Error(TEL_POS_ERR); 120 exit(-1);} 121 else 122 { Telescope_y = (float) atof(value_ptr); 123 Message(TEL_POS_MSG, Telescope_x, Telescope_y);} 124 break; 125 126 case max_events: 127 Message(MAX__EVTS__MSG, max_Events=atol(value_ptr)); 128 break; 129 case energy_cuts: 130 low_Ecut = (float) atof(value_ptr); 131 value_ptr = strtok(NULL, whites); 132 if (value_ptr == NULL) 133 { Error(ENRG_LIMIT_ERR); 134 low_Ecut = 0.; } 135 else 136 { high_Ecut = (float) atof(value_ptr); 137 Message(ENRG_CUTS__MSG, low_Ecut, high_Ecut); } 138 break; 139 case seeds: 140 Seeds[0] = atol(value_ptr); 141 value_ptr = strtok(NULL, whites); 142 if (value_ptr) Seeds[1] = atol(value_ptr); 143 else 144 { Error(SEED_ERROR_ERR); 145 Seeds[0] = 3141592L; } 146 break; 147 case cer_files: 148 ReadCerfiles(parfile); 149 switch_end = TRUE; 150 default: switch_end = TRUE; 151 break; }} 136 152 while (!switch_end); 137 153 138 154 if (filelist == NULL) 139 155 FatalError(FLST_NSPEC_FTL); 140 156 141 157 /* Set random seeds */ … … 143 159 Message(SEED_SET___MSG, Seeds[0], Seeds[1]); 144 160 145 }/* end of ParmsSwitch */161 } /* end of ParmsSwitch */ -
trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h
r725 r1431 6 6 T(output_file), /* output file */ \ 7 7 T(ct_file), /* file with the characteristics of the CT */ \ 8 T(axisdev_file), /* file with the single mirror spot deviations*/ \ 9 T(reflectivity_file), /* file with the mirror reflectivity */ \ 8 10 T(atm_model), /* changes the atmospheric model to be used */ \ 9 11 T(verbose_level), /* defines verbose level of the output */ \ 10 12 T(fixed_target), /* position towards which CT is pointing */ \ 13 T(telescope_position),/* position towards which CT is pointing */ \ 11 14 T(max_events), /* maximum number of event to read */ \ 12 15 T(energy_cuts), /* lowest/highest energy allowed */ \ 13 16 T(seeds), /* seeds for random number generation */ \ 14 T(random_pointing), /* random CT pointing from each shower (hadrons) */ \15 T(repeat_random), /* number of times a random pointing is to be done */ \16 17 T(cer_files) /* start of filename list (must be last) */ 17 18 … … 39 40 #define FIXD_ENABL_MSG /* theta, phi */ \ 40 41 "Using \"fixed_target\" mode with theta=%.2fdeg and phi=%.2fdeg.\n" 42 #define TEL_POS_ERR /* no parms */ \ 43 "Error while parsing \"telescope_position\" fields.\n" \ 44 " *** EXITING! \n" 45 #define TEL_POS_MSG /* theta, phi */ \ 46 "Using \"telescope position\" mode with x=%.1fcm and y=%.1fcm.\n" 41 47 #define MAX__EVTS__MSG /* max. nr. of evts. */ \ 42 48 "Processing at most %ld events.\n" … … 59 65 " *** Please specify some *valid* filename after the cer_files directive.\n" 60 66 67 68 extern char axisdev_filename[256], reflectivity_filename[256]; 69 61 70 #endif -
trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c
r869 r1431 237 237 + 2.0*SQR(rCT[0])*xCT[2] + 2.0*SQR(rCT[1])*xCT[2]); 238 238 239 c = 2*rCT[0]*rCT[2]*x[0]*x[2] + 2*rCT[1]*rCT[2]*x[1]*x[2] 240 - SQR(rCT[2])*SQR(x[0]) - SQR(rCT[2])*SQR(x[1]) 241 - SQR(rCT[0])*SQR(x[2]) - SQR(rCT[1])*SQR(x[2]); 239 /* FIXED Lines below, May 2002, AM : formerly (up to V0.4) 240 * there was a confusion between telescope coordinates xCT and 241 * the original coordinates x. Thanks to T. Hengstebeck for 242 * reporting the bug. 243 */ 244 245 c = 2*rCT[0]*rCT[2]*xCT[0]*xCT[2] + 2*rCT[1]*rCT[2]*xCT[1]*xCT[2] 246 - SQR(rCT[2])*SQR(xCT[0]) - SQR(rCT[2])*SQR(xCT[1]) 247 - SQR(rCT[0])*SQR(xCT[2]) - SQR(rCT[1])*SQR(xCT[2]); 248 242 249 243 250 if ( fabs(a) < 1.e-6 ) { … … 583 590 Debug("@19 xcam-AD %f %f \n", xcam[0], xcam[1]); 584 591 /* CBC */ 585 592 586 593 /* 587 594 ++ -
trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c
r870 r1431 24 24 #include "header.h" 25 25 26 #define MAX(x,y) ((x)>(y)? (x) : (y)) 27 #define MIN(x,y) ((x)>(y)? (y) : (x)) 28 26 29 FILE *rflfile = NULL, /* reflector (output) file */ 27 30 *cerfile = NULL; /* current cerfile */ 28 31 29 32 cphoton *CPhotons = NULL; /* cphoton array */ … … 42 45 static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile); 43 46 44 FILE *chkf = NULL; /**********************/ 45 long myloop; /**********************/ 46 47 void main(void) 47 FILE *chkf = NULL; 48 long myloop; 49 50 Event_end* evt_end; 51 52 main(void) 48 53 { long event = 0L; /* event counter */ 49 54 50 55 /* Read init & geometry parms, init files and vars */ 51 56 init(NULL); 52 chkf = fopen("check", "w"); /**********************/ 57 58 /* chkf = fopen("check", "w"); */ 53 59 54 60 /* Processing loop */ 55 61 while(event < max_Events && GetEvent()) 56 { if (ProcessEvent(cheadp, cerfile, rflfile)) event++; 57 if (event % SHOW_ME == 0) Log(INFO_EVENT_LOG, 58 cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal); } 62 { 63 if (ProcessEvent(cheadp, cerfile, rflfile)) 64 event++; 65 if (event % SHOW_ME == 0) 66 Log(INFO_EVENT_LOG, 67 cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal); 68 } 59 69 60 70 /* Writing final flags */ … … 64 74 Log(RFLF_CLOSE_LOG); 65 75 fclose(rflfile); 66 fclose(chkf); /**********************/76 /* fclose(chkf); */ 67 77 68 78 /* Clean memory and exit */ … … 70 80 Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION)); 71 81 72 }/* end of main */82 } /* end of main */ 73 83 74 84 static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile) … … 82 92 extern void makeOmegaI(float theta, float phi); 83 93 94 char pp[256]; 84 95 85 96 /* Various counters: phs = absphs + refphs[0..3] + cphs */ 86 long phs, /* Number of incoming photons */87 absphs,/* Photons absorbed */88 refphs[4],/* Photons not reflected */89 cphs;/* Number of cphotons */90 91 FILE *tmpf=NULL; /* Temp fp to handle o/f */92 size_t read; /* items read: != 1 on error */93 int overflow=0, /* how many o/f on cphs */94 ref_type,/* ret value from reflection */95 ph;/* photon number */96 97 float first = 1e8f, /* Photon arrival times */98 97 long phs, /* Number of incoming photons */ 98 absphs, /* Photons absorbed */ 99 refphs[4], /* Photons not reflected */ 100 cphs; /* Number of cphotons */ 101 102 FILE *tmpf=NULL; /* Temp fp to handle o/f */ 103 size_t read; /* items read: != 1 on error */ 104 int overflow=0, /* how many o/f on cphs */ 105 ref_type, /* ret value from reflection */ 106 ph; /* photon number */ 107 108 float first = 1e8f, /* Photon arrival times */ 109 last = 0; 99 110 100 111 /* photon quantities */ 101 float wlen, /* photon wavelength */102 theta;/* photon zenith angle */112 float wlen, /* photon wavelength */ 113 theta; /* photon zenith angle */ 103 114 104 115 /* Reset counters and set fileptr */ … … 111 122 /* Calculate OmegaCT matrices */ 112 123 if (is_Fixed_Target) 113 { makeOmega(fixed_Theta, fixed_Phi);124 { makeOmega(fixed_Theta, fixed_Phi); 114 125 makeOmegaI(fixed_Theta, fixed_Phi); } 115 126 else 116 { makeOmega(cheadp->Theta, cheadp->Phi);127 { makeOmega(cheadp->Theta, cheadp->Phi); 117 128 makeOmegaI(cheadp->Theta, cheadp->Phi); } 118 129 memcpy( OmegaCT, Omega, 9*sizeof(float) ); 119 130 memcpy( OmegaICT, OmegaI, 9*sizeof(float) ); 120 131 132 /* ADDED AM May 2002: now we take into account the telescope 133 * position chosen in the Corsika input card via the CERTEL 134 * option, which must be supplied also in the reflector input 135 * card with the option "telescope_position x y". Otherwise 136 * x = y = 0 is assumed, which is the standard mode of 137 * generation for MAGIC. 138 * Here we change rheadp so that the CorePos written to the 139 * .rfl file is referred to the telescope position. However, 140 * I believe that the original "CorePos" vector points 141 * from the core to the origin of coordinates of Corsika, and 142 * therefore after the subtraction of (Telescope_x,Telescope_y) 143 * the resulting vector CorePos points from the core to the 144 * telescope! 145 */ 146 147 rheadp->CorePos[0][0] += Telescope_x; 148 rheadp->CorePos[1][0] += Telescope_y; 149 121 150 /* Loop on data blocks */ 122 while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile)) 123 && strncmp((char *)Photons, "EVTE", 4)) 124 125 /* Loop on phs inside block: exit when wlength is 0 */ 126 { for (ph=0; ph<PH_IN_DATABLOCK; ph++) 127 { if (Photons[ph].w == 0) break; 151 152 while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile))) 153 /* Loop on phs inside block: exit when wlength is 0 */ 154 { 155 /* If "event end" flag is found, read relevant quantities 156 * from event end subblock (added June 2002, AM): 157 */ 158 159 if (strncmp((char *)Photons, "EVTE", 4) == 0 ) 160 { 161 evt_end = (Event_end*) Photons; 162 rheadp->longi_Nmax = evt_end->longi_Nmax; 163 rheadp->longi_t0 = evt_end->longi_t0; 164 rheadp->longi_tmax = evt_end->longi_tmax; 165 rheadp->longi_a = evt_end->longi_a; 166 rheadp->longi_b = evt_end->longi_b; 167 rheadp->longi_c = evt_end->longi_c; 168 rheadp->longi_chi2 = evt_end->longi_chi2; 169 break; 170 } 171 172 for (ph=0; ph<PH_IN_DATABLOCK; ph++) 173 { 174 if (Photons[ph].w <= 0.) break; 128 175 129 176 CPhotons[cphs].w = Photons[ph].w; 130 177 Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.); 131 Photons[ph].x -= cheadp->CorePos[0][0]; 132 Photons[ph].y -= cheadp->CorePos[1][0]; 178 179 /* TEMPORARY FIX, AM Nov 2001: we found that sometimes the value 180 stored in Photons[ph].w is not correct, and can result in a 181 wavelength beyond 600 nm, which makes the program crash later. 182 Now we force wlen to its expected range: 183 */ 184 185 wlen = MIN(MAX(290.,wlen),600.); 186 187 188 /* ADDED AM May 2002: now we take into account the telescope 189 * position chosen in the Corsika input card via the CERTEL 190 * option, which must be supplied also in the reflector input 191 * card with the option "telescope_position x y". Otherwise 192 * x = y = 0 is assumed, which is the standard mode of 193 * generation for MAGIC. 194 */ 195 196 Photons[ph].x -= cheadp->CorePos[0][0]+Telescope_x; 197 Photons[ph].y -= cheadp->CorePos[1][0]+Telescope_y; 133 198 134 199 /* Increment number of photons */ … … 137 202 /* Calculate some quantities */ 138 203 theta = (float) acos(sqrt( 139 140 - Photons[ph].v*Photons[ph].v));204 MAX(0., 1.f - Photons[ph].u*Photons[ph].u 205 - Photons[ph].v*Photons[ph].v))); 141 206 142 207 /* Check absorption */ 208 209 143 210 if (absorption(wlen, Photons[ph].h, theta)) 144 211 absphs++; 145 212 146 213 /* Check reflection */ 147 214 else if (0 != (ref_type = 148 215 ph2cph(&Photons[ph], &CPhotons[cphs]))) 149 150 else 151 {152 153 Photons[ph].x,Photons[ph].y,154 Photons[ph].u,Photons[ph].v,theta);155 156 CPhotons[cphs].y);216 refphs[ref_type-1]++; 217 else /* Photon passed */ 218 { 219 Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph, 220 Photons[ph].x,Photons[ph].y, 221 Photons[ph].u,Photons[ph].v,theta); 222 Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x, 223 CPhotons[cphs].y); 157 224 158 225 … … 162 229 163 230 /* Update cphs */ 164 if (++cphs == NR_OF_CPHOTONS) 165 {231 if (++cphs == NR_OF_CPHOTONS) /* Overflow */ 232 { 166 233 /* if it is the first o/f open a tempfile */ 167 234 if (overflow++ == 0) 168 { Log("Spooling... ");235 { Log("Spooling... "); 169 236 tmpf = tmpfile(); 170 237 if (tmpf == NULL) FatalError(TEMP_ERROR_FTL); } … … 176 243 cphs = 0; 177 244 178 }/* if overflow */179 }/* else (=Photon passed) */180 }/* end of loop inside datablock */181 }/* end of loop on datablocks */245 } /* if overflow */ 246 } /* else (=Photon passed) */ 247 } /* end of loop inside datablock */ 248 } /* end of loop on datablocks */ 182 249 183 250 /* Check if there was an error while reading cerfile */ 184 251 if (read != 1) fseek(rflfile, writep, SEEK_SET); 185 252 186 else 187 188 {/* Write "start of event" flag */253 else /* no error: write the new event */ 254 255 { /* Write "start of event" flag */ 189 256 fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile); 190 257 … … 202 269 rheadp->CPhotons = (long) overflow * NR_OF_CPHOTONS + cphs; 203 270 fwrite(rheadp, sizeof(RflHeader), 1, rflfile); 204 for (myloop=0; myloop<sizeof(RflHeader)/4; myloop++) 205 fprintf(chkf, "%e ", *((float *)rheadp+myloop)); 206 fputc('\n', chkf); 271 272 /* for (myloop=0; myloop<sizeof(RflHeader)/4; myloop++) 273 * fprintf(chkf, "%e ", *((float *)rheadp+myloop)); 274 * fputc('\n', chkf); 275 */ 207 276 208 277 /* If there was an overflow, append data from tempfile */ 209 278 if (overflow) 210 {279 { 211 280 /* Unload data from CPhotons */ 212 281 fwrite(CPhotons, sizeof(cphoton), cphs, tmpf); 213 282 214 283 /* Transfer data */ 215 fseek(tmpf, 0L, SEEK_SET); 284 fseek(tmpf, 0L, SEEK_SET); /* Start from the beginning */ 216 285 while (overflow--) 217 { fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);286 { fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf); 218 287 fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile); } 219 288 … … 236 305 237 306 return read == 1; 238 }/* end of ProcessEvent */307 } /* end of ProcessEvent */ 239 308 240 309 static int GetEvent(void) 241 310 { int found = FALSE, /* event found */ 242 isWrong = FALSE;/* cerfile is wrong */243 static int newFile = TRUE;/* if TRUE, check if cerfile is valid */244 245 do246 { /*In case the just-opened file is a valid cerfile,247 starting with a RUNH, loop until a valid event is found:248 id est with: first_Event <= EvtNumber <= last_Event249 and low_Ecut <= Etotal <= high_Ecut250 If the search was successful, "found" is set to TRUE.251 If there are reading errors, "isWrong" is set to TRUE. */311 isWrong = FALSE; /* cerfile is wrong */ 312 static int newFile = TRUE; /* if TRUE, check if cerfile is valid */ 313 314 do 315 { /* In case the just-opened file is a valid cerfile, 316 starting with a RUNH, loop until a valid event is found: 317 id est with: first_Event <= EvtNumber <= last_Event 318 and low_Ecut <= Etotal <= high_Ecut 319 If the search was successful, "found" is set to TRUE. 320 If there are reading errors, "isWrong" is set to TRUE. */ 252 321 253 322 if (newFile … … 308 377 309 378 FILE *GetNextFile(char *cername) 310 { FILE * f= NULL; /* return value (cerfile ptr) */379 { FILE *inputfile = NULL; /* return value (cerfile ptr) */ 311 380 char *value_ptr; /* ptr at parm value */ 312 381 extern char line[]; /* white chars (init) */ … … 324 393 325 394 /* If you found a line with some meaning, try to open the file */ 326 if (( f=fopen(value_ptr, "r")) == NULL)395 if ((inputfile=fopen(value_ptr, "r")) == NULL) 327 396 Message(CERF_NFND__MSG, value_ptr); 328 397 … … 346 415 last_Event = 1000000; }}} 347 416 348 while ( f== NULL); /* Loop until a readable file is found */417 while (inputfile == NULL); /* Loop until a readable file is found */ 349 418 350 419 /* If a new cerfile is met, write "START_OF_RUN" */ 351 if ( f) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);352 353 return f;420 if (inputfile) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile); 421 422 return inputfile; 354 423 } /* end of GetNextFile */ 355 424 … … 378 447 { free(CPhotons); Log(VECT_FREE__LOG, "CPhotons"); } 379 448 } /* end of clean */ 449 -
trunk/MagicSoft/Simulation/Detector/ReflectorII/version.h
r725 r1431 6 6 7 7 #define PROGRAM reflector 8 #define VERSION 0. 48 #define VERSION 0.5 9 9 10 10 #endif
Note:
See TracChangeset
for help on using the changeset viewer.
