1 | /* ======================================================================== *\
|
---|
2 | !
|
---|
3 | ! *
|
---|
4 | ! * This file is part of CheObs, the Modular Analysis and Reconstruction
|
---|
5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
8 | ! *
|
---|
9 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
10 | ! * documentation for any purpose is hereby granted without fee,
|
---|
11 | ! * provided that the above copyright notice appears in all copies and
|
---|
12 | ! * that both that copyright notice and this permission notice appear
|
---|
13 | ! * in supporting documentation. It is provided "as is" without express
|
---|
14 | ! * or implied warranty.
|
---|
15 | ! *
|
---|
16 | !
|
---|
17 | !
|
---|
18 | ! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
19 | !
|
---|
20 | ! Copyright: CheObs Software Development, 2000-2009
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | //////////////////////////////////////////////////////////////////////////////
|
---|
26 | //
|
---|
27 | // MSimAtmosphere
|
---|
28 | //
|
---|
29 | // Task to calculate wavelength or incident angle dependent absorption
|
---|
30 | //
|
---|
31 | // Input Containers:
|
---|
32 | // MPhotonEvent
|
---|
33 | // MCorsikaRunHeader
|
---|
34 | //
|
---|
35 | // Output Containers:
|
---|
36 | // MPhotonEvent
|
---|
37 | //
|
---|
38 | //////////////////////////////////////////////////////////////////////////////
|
---|
39 | #include "MSimAtmosphere.h"
|
---|
40 |
|
---|
41 | #include <fstream>
|
---|
42 |
|
---|
43 | #include <TGraph.h>
|
---|
44 | #include <TRandom.h>
|
---|
45 |
|
---|
46 | #include "MLog.h"
|
---|
47 | #include "MLogManip.h"
|
---|
48 |
|
---|
49 | #include "MParList.h"
|
---|
50 |
|
---|
51 | #include "MCorsikaRunHeader.h"
|
---|
52 | #include "MPhotonEvent.h"
|
---|
53 | #include "MPhotonData.h"
|
---|
54 |
|
---|
55 | ClassImp(MSimAtmosphere);
|
---|
56 |
|
---|
57 | using namespace std;
|
---|
58 |
|
---|
59 | // ==========================================================================
|
---|
60 | //
|
---|
61 | // January 2002, A. Moralejo: We now precalculate the slant paths for the
|
---|
62 | // aerosol and Ozone vertical profiles, and then do an interpolation in
|
---|
63 | // wavelength for every photon to get the optical depths. The parameters
|
---|
64 | // used, defined below, have been taken from "Atmospheric Optics", by
|
---|
65 | // L. Elterman and R.B. Toolin, chapter 7 of the "Handbook of geophysics
|
---|
66 | // and Space environments". (S.L. Valley, editor). McGraw-Hill, NY 1965.
|
---|
67 | //
|
---|
68 | // WARNING: the Mie scattering and the Ozone absorption are implemented
|
---|
69 | // to work only with photons produced at a height a.s.l larger than the
|
---|
70 | // observation level. So this is not expected to work well for simulating
|
---|
71 | // the telescope pointing at theta > 90 deg (for instance for neutrino
|
---|
72 | // studies. Rayleigh scattering works even for light coming from below.
|
---|
73 | //
|
---|
74 | // Fixed bugs (of small influence) in Mie absorption implementation: there
|
---|
75 | // were errors in the optical depths table, as well as a confusion:
|
---|
76 | // height a.s.l. was used as if it was height above the telescope level.
|
---|
77 | // The latter error was also present in the Ozone aborption implementation.
|
---|
78 | //
|
---|
79 | // On the other hand, now we have the tables AE_ABI and OZ_ABI with optical
|
---|
80 | // depths between sea level and a height h (before it was between 2km a.s.l
|
---|
81 | // and a height h). So that we can simulate also in the future a different
|
---|
82 | // observation level.
|
---|
83 | //
|
---|
84 | // AM: WARNING: IF VERY LARGE zenith angle simulations are to be done (say
|
---|
85 | // above 85 degrees, for neutrino primaries or any other purpose) this code
|
---|
86 | // will have to be adapted accordingly and checked, since in principle it has
|
---|
87 | // been written and tested to simulate the absorption of Cherenkov photons
|
---|
88 | // arriving at the telescope from above.
|
---|
89 | //
|
---|
90 | // AM: WARNING 2: not to be used for wavelengths outside the range 250-700 nm
|
---|
91 | //
|
---|
92 | // January 2003, Abelardo Moralejo: found error in Ozone absorption treatment.
|
---|
93 | // At large zenith angles, the air mass used was the one calculated for
|
---|
94 | // Rayleigh scattering, but since the Ozone distribution is rather different
|
---|
95 | // from the global distribution of air molecules, this is not a good
|
---|
96 | // approximation. Now I have left in this code only the Rayleigh scattering,
|
---|
97 | // and moved to atm.c the Mie scattering and Ozone absorption calculated in
|
---|
98 | // a better way.
|
---|
99 | //
|
---|
100 | // A. Moralejo, January 2003: added some parameters for Mie scattering
|
---|
101 | // and Ozone absorption derived from the clear standard atmosphere model
|
---|
102 | // in "Atmospheric Optics", by L. Elterman and R.B. Toolin, chapter 7 of
|
---|
103 | // the "Handbook of geophysics and Space environments". S.L. Valley,
|
---|
104 | // editor. McGraw-Hill, NY 1965.
|
---|
105 | //
|
---|
106 | // AM, Jan 2003: Changed the meaning of the argument height: now it is the
|
---|
107 | // true height above sea level at which a photon has been emitted, before
|
---|
108 | // it was the height given by Corsika, measured in the vertical of the
|
---|
109 | // observer and not in the vertical of the emitting particle.
|
---|
110 | //
|
---|
111 | //
|
---|
112 | // MAGIC-Winter and MAGIC-Summer by M. Haffke,
|
---|
113 | // parametrizing profiles obtained with MSIS:
|
---|
114 | // http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm
|
---|
115 | //
|
---|
116 | //
|
---|
117 | // The MAGIC-Winter and MAGIC-Summer parametrisations reproduce the MSIS
|
---|
118 | // profiles for the 3 atmospheric layers from 0 up to 40 km height. Beyond
|
---|
119 | // that height, it was not possible to achieve a good fit, but the amount
|
---|
120 | // of residual atmosphere is so small that light absorption would be
|
---|
121 | // negligible anyway. Showers develop well below 40 km.
|
---|
122 | //
|
---|
123 | //
|
---|
124 | // The mass overburden is given by T = AATM + BATM * exp(-h/CATM)
|
---|
125 | // The last layer of the US standard atmosphere (in which T varies
|
---|
126 | // linearly with h) is above 100 km and has not been included here
|
---|
127 | // because it should not matter.
|
---|
128 | //
|
---|
129 | class MAtmRayleigh
|
---|
130 | {
|
---|
131 | private:
|
---|
132 | static const Double_t fgMeanFreePath; // [g/cm^2] Mean free path for scattering Rayleigh XR
|
---|
133 |
|
---|
134 | Double_t fR; // [cm] Earth radius to be used
|
---|
135 |
|
---|
136 | Double_t fHeight[5]; // Layer boundaries (atmospheric layer)
|
---|
137 |
|
---|
138 | // Parameters of the different atmospheres. We use the same parametrization
|
---|
139 | // shape as in Corsika atmospheric models (see Corsika manual, appendix D).
|
---|
140 | // The values here can be/are obtained from the Corsika output
|
---|
141 | //Float_t fAtmA[4]; // The index refers to the atmospheric layer (starting from sea level and going upwards)
|
---|
142 | Float_t fAtmB[4]; // The index refers to the atmospheric layer (starting from sea level and going upwards)
|
---|
143 | Float_t fAtmC[4]; // The index refers to the atmospheric layer (starting from sea level and going upwards)
|
---|
144 |
|
---|
145 | Double_t fRho[5]; // Precalculated integrals for rayleigh scatterning
|
---|
146 |
|
---|
147 | // --------------------------------------------------------------------------
|
---|
148 | //
|
---|
149 | // Precalcalculate the integrals from the observer level to the next
|
---|
150 | // atmpsheric layer below including the lower boundary. Thus a
|
---|
151 | // correct calculation is reduced to the calculation of the upper
|
---|
152 | // boundary.
|
---|
153 | //
|
---|
154 | // fRho[0] = B0;
|
---|
155 | // fRho[1] = B0-A0 + B1;
|
---|
156 | // fRho[2] = B0-A0 + B1-A1 + B2;
|
---|
157 | // fRho[3] = B0-A0 + B1-A1 + B2+A2 + B3;
|
---|
158 | // fRho[4] = B0-A0 + B1-A1 + B2+A2 + B3 - A3;
|
---|
159 | //
|
---|
160 | void PreCalcRho()
|
---|
161 | {
|
---|
162 | // Limits (height in cm) of the four layers in which
|
---|
163 | // atmosphere is parametrized.
|
---|
164 | // This is a stupid trick giving 0 for the integrals below
|
---|
165 | // the observer
|
---|
166 |
|
---|
167 | const Double_t h[5] =
|
---|
168 | {
|
---|
169 | fObsLevel, // 0km
|
---|
170 | TMath::Max(fObsLevel, 4e5), // 4km
|
---|
171 | 1.0e6, // 10km
|
---|
172 | 4.0e6, // 40km
|
---|
173 | 1.0e7 // 100km
|
---|
174 | };
|
---|
175 |
|
---|
176 | memcpy(fHeight, h, sizeof(Double_t)*5);
|
---|
177 |
|
---|
178 | fRho[0] = 0;
|
---|
179 | for (int i=0; i<4; i++)
|
---|
180 | {
|
---|
181 | const Double_t b = fAtmB[i];
|
---|
182 | const Double_t c = fAtmC[i];
|
---|
183 |
|
---|
184 | const Double_t h1 = h[i+1];
|
---|
185 | const Double_t h0 = h[i];
|
---|
186 |
|
---|
187 | const Double_t B = b*TMath::Exp(-h0/c);
|
---|
188 | const Double_t A = b*TMath::Exp(-h1/c);
|
---|
189 |
|
---|
190 | // Calculate rho for the i-th layer from the lower
|
---|
191 | // to the higher layer boundary.
|
---|
192 | // If height is within the layer only calculate up to height.
|
---|
193 | fRho[i] += B;
|
---|
194 | fRho[i+1] = fRho[i] - A;
|
---|
195 | }
|
---|
196 | }
|
---|
197 |
|
---|
198 | protected:
|
---|
199 | Double_t fObsLevel; // [cm] observation level a.s.l.
|
---|
200 |
|
---|
201 | public:
|
---|
202 | // Init an atmosphere from the data stored in MCorsikaRunHeader
|
---|
203 | MAtmRayleigh(const MCorsikaRunHeader &h)
|
---|
204 | {
|
---|
205 | Init(h);
|
---|
206 | }
|
---|
207 |
|
---|
208 | // Defualt constructor
|
---|
209 | MAtmRayleigh() : fObsLevel(-1) { }
|
---|
210 |
|
---|
211 | // Check if the ovservation level has been correctly initialized
|
---|
212 | // Used as a marker for correct initialization
|
---|
213 | Bool_t IsValid() const { return fObsLevel>=0; }
|
---|
214 |
|
---|
215 | // Get the Earth radius to be used
|
---|
216 | Double_t R() const { return fR; }
|
---|
217 |
|
---|
218 | // Init an atmosphere from the data stored in MCorsikaRunHeader
|
---|
219 | // This initialized fObsLevel, fR, fAtmB and fAtmC and
|
---|
220 | // PreCalcRho
|
---|
221 | void Init(const MCorsikaRunHeader &h)
|
---|
222 | {
|
---|
223 | // Observation level above earth radius
|
---|
224 | fObsLevel = h.GetObsLevel();
|
---|
225 |
|
---|
226 | // Use earth radius as defined in Corsika
|
---|
227 | fR = h.EarthRadius();
|
---|
228 |
|
---|
229 | //memcpy(fAtmA, (Float_t*)h.GetAtmosphericCoeffA(), sizeof(Float_t)*4);
|
---|
230 | memcpy(fAtmB, (Float_t*)h.GetAtmosphericCoeffB(), sizeof(Float_t)*4);
|
---|
231 | memcpy(fAtmC, (Float_t*)h.GetAtmosphericCoeffC(), sizeof(Float_t)*4);
|
---|
232 |
|
---|
233 | PreCalcRho();
|
---|
234 | }
|
---|
235 |
|
---|
236 | // Return the vertical thickness between the observer and height.
|
---|
237 | // Therefor the integral of the layers below (including the lower
|
---|
238 | // boudary) Have been precalculated by PreCalcRho
|
---|
239 | Double_t GetVerticalThickness(Double_t height) const
|
---|
240 | {
|
---|
241 | // FIXME: We could store the start point above obs-level
|
---|
242 | // (Does this really gain anything?)
|
---|
243 | Int_t i=0;
|
---|
244 | while (i<4 && height>fHeight[i+1])
|
---|
245 | i++;
|
---|
246 |
|
---|
247 | const Double_t b = fAtmB[i];
|
---|
248 | const Double_t c = fAtmC[i];
|
---|
249 |
|
---|
250 | return fRho[i] - b*TMath::Exp(-height/c);
|
---|
251 | }
|
---|
252 |
|
---|
253 | /*
|
---|
254 | // The "orginal" code for the vertical thickness
|
---|
255 | Double_t GetVerticalThickness(Double_t obslev, Double_t height) const
|
---|
256 | {
|
---|
257 | // This is a C++-version of the original code from attenu.c
|
---|
258 |
|
---|
259 | // Limits (height in cm) of the four layers in which atmosphere is parametrized:
|
---|
260 | const double lahg[5] =
|
---|
261 | {
|
---|
262 | obslev,
|
---|
263 | TMath::Max(obslev, 4e5),
|
---|
264 | 1.0e6,
|
---|
265 | 4.0e6,
|
---|
266 | 1.0e7
|
---|
267 | };
|
---|
268 |
|
---|
269 | Double_t Rho_Tot = 0.0;
|
---|
270 | for (int i=0; i<4; i++)
|
---|
271 | {
|
---|
272 | const Double_t b = fAtmB[i];
|
---|
273 | const Double_t c = fAtmC[i];
|
---|
274 |
|
---|
275 | const Double_t h0 = TMath::Min(height, lahg[i+1]);
|
---|
276 | const Double_t h1 = lahg[i];
|
---|
277 |
|
---|
278 | // Calculate rho for the i-th layer from the lower
|
---|
279 | // to the higher layer boundary.
|
---|
280 | // If height is within the layer only calculate up to height.
|
---|
281 | Rho_Tot += b*(exp(-h1/c) - exp(-h0/c));
|
---|
282 |
|
---|
283 | if (lahg[i+1] > height)
|
---|
284 | break;
|
---|
285 | }
|
---|
286 |
|
---|
287 | return Rho_Tot;
|
---|
288 | }
|
---|
289 | */
|
---|
290 | Double_t CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const
|
---|
291 | {
|
---|
292 | // sin2: sin(theta)^2
|
---|
293 | // height is the true height a.s.l.
|
---|
294 |
|
---|
295 | // LARGE ZENITH ANGLE FACTOR (AIR MASS FACTOR):
|
---|
296 | // Air mass factor "airmass" calculated using a one-exponential
|
---|
297 | // density profile for the atmosphere,
|
---|
298 | //
|
---|
299 | // rho = rho_0 exp(-height/hscale) with hscale = 7.4 km
|
---|
300 | //
|
---|
301 | // The air mass factor is defined as I(theta)/I(0), the ratio of
|
---|
302 | // the optical paths I (in g/cm2) traversed between two given
|
---|
303 | // heights, at theta and at 0 deg z.a.
|
---|
304 |
|
---|
305 | const Double_t H = height-fObsLevel;
|
---|
306 | const Double_t h = 2*H;
|
---|
307 |
|
---|
308 | // Scale-height (cm) for Exponential density profile
|
---|
309 | const Double_t hscale = 7.4e5;
|
---|
310 | const Double_t f = 2*hscale;
|
---|
311 |
|
---|
312 | // Precalc R*cos(theta)^2 (FIXME: Is ph.GetCosW2 more precise?)
|
---|
313 | const Double_t Rcos2 = fR * (1-sin2); // cos2 = 1 - sin2
|
---|
314 |
|
---|
315 | const Double_t x1 = TMath::Sqrt((Rcos2 ) / f);
|
---|
316 | const Double_t x2 = TMath::Sqrt((Rcos2 + 2*h) / f);
|
---|
317 | const Double_t x3 = TMath::Sqrt((fR ) / f);
|
---|
318 | const Double_t x4 = TMath::Sqrt((fR + 2*h) / f);
|
---|
319 |
|
---|
320 | // Return a -1 transmittance in the case the photon comes
|
---|
321 | // exactly from the observation level, because in that case the
|
---|
322 | // "air mass factor" would become infinity and the calculation
|
---|
323 | // is not valid!
|
---|
324 | if (fabs(x3-x4) < 1.e-10)
|
---|
325 | return -1.;
|
---|
326 |
|
---|
327 | const Double_t e12 = erfc(x1) - erfc(x2);
|
---|
328 | const Double_t e34 = erfc(x3) - erfc(x4);
|
---|
329 |
|
---|
330 | const Double_t airmass = TMath::Exp(-fR*sin2 / f) * e12/e34;
|
---|
331 |
|
---|
332 | // Calculate the traversed "vertical thickness" of air using the
|
---|
333 | // US Standard atmosphere:
|
---|
334 | const Double_t Rho_tot = GetVerticalThickness(/*fObsLevel,*/ height);
|
---|
335 |
|
---|
336 | // We now convert from "vertical thickness" to "slanted thickness"
|
---|
337 | // traversed by the photon on its way to the telescope, simply
|
---|
338 | // multiplying by the air mass factor m:
|
---|
339 | const Double_t Rho_Fi = airmass * Rho_tot;
|
---|
340 |
|
---|
341 | // Finally we calculate the transmission coefficient for the Rayleigh
|
---|
342 | // scattering:
|
---|
343 | // AM Dec 2002, introduced ABS below to account (in the future) for
|
---|
344 | // possible photons coming from below the observation level.
|
---|
345 |
|
---|
346 | const Double_t wl = 400./wavelength;
|
---|
347 |
|
---|
348 | // Mean free path for scattering Rayleigh XR (g/cm^2)
|
---|
349 | return TMath::Exp(-TMath::Abs(Rho_Fi/fgMeanFreePath)*wl*wl*wl*wl);
|
---|
350 | }
|
---|
351 | };
|
---|
352 |
|
---|
353 | // ==========================================================================
|
---|
354 |
|
---|
355 | class MAtmosphere : public MAtmRayleigh
|
---|
356 | {
|
---|
357 | private:
|
---|
358 | static const Double_t STEPTHETA; // aprox. 1 degree
|
---|
359 |
|
---|
360 | // Aerosol number density for 31 heights a.s.l., from 0 to 30 km,
|
---|
361 | // in 1 km steps (units: cm^-3)
|
---|
362 | static const Double_t aero_n[31];
|
---|
363 |
|
---|
364 | // Ozone concentration for 51 heights a.s.l., from 0 to 50 km,
|
---|
365 | // in 1 km steps (units: cm/km)
|
---|
366 | static const Double_t oz_conc[51];
|
---|
367 |
|
---|
368 | // aerosol_path contains the path integrals for the aerosol number
|
---|
369 | // density (relative to the number density at sea level) between the
|
---|
370 | // observation level and a height h for different zenith angles. The
|
---|
371 | // first index indicate height above sea level in units of 100m, the
|
---|
372 | // second is the zenith angle in degrees.
|
---|
373 | float aerosol_path[301][90];
|
---|
374 |
|
---|
375 | // ozone_path contains the path integrals for the ozone concentration
|
---|
376 | // between the observation level and a height h for different zenith
|
---|
377 | // angles. The first index indicate height above sea level in units
|
---|
378 | // of 100m, the second is the zenith angle in degrees.
|
---|
379 | float ozone_path[501][90];
|
---|
380 |
|
---|
381 | // Interpolate the graph at wavelength
|
---|
382 | Double_t GetBeta(Double_t wavelength, const TGraph &g) const
|
---|
383 | {
|
---|
384 | // FIXME: This might not be the fastest because range
|
---|
385 | // checks are done for each call!
|
---|
386 | return g.GetN()==0 ? 0 : g.Eval(wavelength)*1e-5; // from km^-1 to cm^-1
|
---|
387 | /*
|
---|
388 | // Linear interpolation
|
---|
389 | // (FIXME: Is it faster to be replaced with a binary search?)
|
---|
390 | // ( This might be faster because we have more photons
|
---|
391 | // with smaller wavelengths)
|
---|
392 | //int index;
|
---|
393 | //for (index = 1; index <g.GetN()-1; index++)
|
---|
394 | // if (wavelength < g.GetX()[index])
|
---|
395 | // break;
|
---|
396 | const Int_t index = TMath::BinarySearch(g.GetN(), g.GetX(), wavelength)+1;
|
---|
397 |
|
---|
398 | const Double_t t0 = g.GetY()[index-1];
|
---|
399 | const Double_t t1 = g.GetY()[index];
|
---|
400 |
|
---|
401 | const Double_t w0 = g.GetX()[index-1];
|
---|
402 | const Double_t w1 = g.GetX()[index];
|
---|
403 |
|
---|
404 | const Double_t beta0 = t0+(t1-t0)*(wavelength-w0)/(w1-w0);
|
---|
405 |
|
---|
406 | return beta0 * 1e-5; // from km^-1 to cm^-1
|
---|
407 | */
|
---|
408 | }
|
---|
409 |
|
---|
410 |
|
---|
411 | //MSpline3 *fAbsCoeffOzone;
|
---|
412 | //MSpline3 *fAbsCoeffAerosols;
|
---|
413 |
|
---|
414 | TGraph *fAbsCoeffOzone;
|
---|
415 | TGraph *fAbsCoeffAerosols;
|
---|
416 |
|
---|
417 | public:
|
---|
418 | MAtmosphere(const MCorsikaRunHeader &h) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
|
---|
419 | {
|
---|
420 | Init(h, "ozone.txt", "aerosols.txt");
|
---|
421 | }
|
---|
422 |
|
---|
423 | MAtmosphere(const char *name1="ozone.txt", const char *name2="aerosols.txt") : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
|
---|
424 | {
|
---|
425 | InitOzone(name1);
|
---|
426 | InitAerosols(name2);
|
---|
427 | }
|
---|
428 |
|
---|
429 | ~MAtmosphere()
|
---|
430 | {
|
---|
431 | if (fAbsCoeffOzone)
|
---|
432 | delete fAbsCoeffOzone;
|
---|
433 | if (fAbsCoeffAerosols)
|
---|
434 | delete fAbsCoeffAerosols;
|
---|
435 | }
|
---|
436 |
|
---|
437 | Float_t GetWavelengthMin() const { return fAbsCoeffOzone && fAbsCoeffAerosols ? TMath::Max(fAbsCoeffOzone->GetX()[0], fAbsCoeffAerosols->GetX()[0]) : -1; }
|
---|
438 | Float_t GetWavelengthMax() const { return fAbsCoeffOzone && fAbsCoeffAerosols ? TMath::Min(fAbsCoeffOzone->GetX()[fAbsCoeffOzone->GetN()-1], fAbsCoeffAerosols->GetX()[fAbsCoeffAerosols->GetN()-1]) : -1; }
|
---|
439 |
|
---|
440 | Bool_t HasValidOzone() const { return fAbsCoeffOzone && fAbsCoeffOzone->GetN()>0; }
|
---|
441 | Bool_t HasValidAerosol() const { return fAbsCoeffAerosols && fAbsCoeffAerosols->GetN()>0; }
|
---|
442 |
|
---|
443 | Bool_t IsAllValid() const { return IsValid() && HasValidOzone() && HasValidAerosol(); }
|
---|
444 |
|
---|
445 | void PreCalcOzone()
|
---|
446 | {
|
---|
447 | // It follows a precalculation of the slant path integrals we need
|
---|
448 | // for the estimate of the Mie scattering and Ozone absorption:
|
---|
449 | Double_t dh = 1.e3;
|
---|
450 | const Double_t STEPTHETA = 1.74533e-2; // aprox. 1 degree
|
---|
451 |
|
---|
452 | // Ozone absorption
|
---|
453 | for (Int_t j = 0; j < 90; j++)
|
---|
454 | {
|
---|
455 | const Double_t theta = j * STEPTHETA;
|
---|
456 | const Double_t sin2 = sin(theta)*sin(theta);
|
---|
457 | const Double_t H = R()+fObsLevel;
|
---|
458 |
|
---|
459 | Double_t path_slant = 0;
|
---|
460 | for (Double_t h = fObsLevel; h <= 50e5; h += dh)
|
---|
461 | {
|
---|
462 | // h is the true height vertical above ground
|
---|
463 | if (fmod(h,1e4) == 0)
|
---|
464 | ozone_path[(int)(h/1e4)][j] = path_slant;
|
---|
465 |
|
---|
466 | const Double_t km = h/1e5;
|
---|
467 | const Int_t i = TMath::FloorNint(km);
|
---|
468 | const Double_t l = R()+h;
|
---|
469 |
|
---|
470 | const Double_t L = TMath::Sqrt(l*l - H*H * sin2);
|
---|
471 | const Double_t f = dh * l / L;
|
---|
472 |
|
---|
473 | // Linear interpolation at h/1e5
|
---|
474 | Double_t interpol = oz_conc[i] + fmod(km, 1) * (oz_conc[i+1]-oz_conc[i]);
|
---|
475 |
|
---|
476 | path_slant += f * interpol;
|
---|
477 | }
|
---|
478 | }
|
---|
479 | }
|
---|
480 |
|
---|
481 | void PreCalcAerosol()
|
---|
482 | {
|
---|
483 | // It follows a precalculation of the slant path integrals we need
|
---|
484 | // for the estimate of the Mie scattering and Ozone absorption:
|
---|
485 | Double_t dh = 1.e3;
|
---|
486 | const Double_t STEPTHETA = 1.74533e-2; // aprox. 1 degree
|
---|
487 |
|
---|
488 | /* Mie (aerosol): */
|
---|
489 | for (Int_t j = 0; j < 90; j++)
|
---|
490 | {
|
---|
491 | const Double_t theta = j * STEPTHETA;
|
---|
492 | const Double_t sin2 = sin(theta)*sin(theta);
|
---|
493 | const Double_t H = R()+fObsLevel;
|
---|
494 |
|
---|
495 | Double_t path_slant = 0;
|
---|
496 | for (Double_t h = fObsLevel; h <= 30e5; h += dh)
|
---|
497 | {
|
---|
498 | // h is the true height vertical above ground
|
---|
499 | if (fmod(h,1e4) == 0)
|
---|
500 | aerosol_path[(int)(h/1e4)][j] = path_slant;
|
---|
501 |
|
---|
502 | const Double_t km = h/1e5;
|
---|
503 | const Int_t i = TMath::FloorNint(km);
|
---|
504 | const Double_t l = R()+h;
|
---|
505 |
|
---|
506 | const Double_t L = TMath::Sqrt(l*l - H*H * sin2);
|
---|
507 | const Double_t f = dh * l / L;
|
---|
508 |
|
---|
509 | // Linear interpolation at h/1e5
|
---|
510 | Double_t interpol = aero_n[i] + fmod(km, 1)*(aero_n[i+1]-aero_n[i]);
|
---|
511 |
|
---|
512 | path_slant += f * interpol/aero_n[0]; // aero_n [km^-1]
|
---|
513 | }
|
---|
514 | }
|
---|
515 | }
|
---|
516 |
|
---|
517 | Bool_t InitOzone(const TString name="")
|
---|
518 | {
|
---|
519 | if (!name.IsNull())
|
---|
520 | {
|
---|
521 | if (fAbsCoeffOzone)
|
---|
522 | delete fAbsCoeffOzone;
|
---|
523 |
|
---|
524 | fAbsCoeffOzone = new TGraph(name);
|
---|
525 | fAbsCoeffOzone->Sort();
|
---|
526 | }
|
---|
527 |
|
---|
528 | if (!HasValidAerosol())
|
---|
529 | return kFALSE;
|
---|
530 |
|
---|
531 | if (IsValid())
|
---|
532 | PreCalcOzone();
|
---|
533 |
|
---|
534 | return kTRUE;
|
---|
535 | }
|
---|
536 |
|
---|
537 | Bool_t InitAerosols(const TString name="")
|
---|
538 | {
|
---|
539 | if (!name.IsNull())
|
---|
540 | {
|
---|
541 | if (fAbsCoeffAerosols)
|
---|
542 | delete fAbsCoeffAerosols;
|
---|
543 |
|
---|
544 | fAbsCoeffAerosols = new TGraph(name);
|
---|
545 | fAbsCoeffAerosols->Sort();
|
---|
546 | }
|
---|
547 |
|
---|
548 | if (!HasValidAerosol())
|
---|
549 | return kFALSE;
|
---|
550 |
|
---|
551 | if (IsValid())
|
---|
552 | PreCalcOzone();
|
---|
553 |
|
---|
554 | return kTRUE;
|
---|
555 | }
|
---|
556 |
|
---|
557 | void Init(const MCorsikaRunHeader &h, const char *name1=0, const char *name2=0)
|
---|
558 | {
|
---|
559 | MAtmRayleigh::Init(h);
|
---|
560 |
|
---|
561 | InitOzone(name1);
|
---|
562 | InitAerosols(name2);
|
---|
563 | }
|
---|
564 | /*
|
---|
565 | Double_t GetOz(Double_t height, Double_t theta) const
|
---|
566 | {
|
---|
567 | // Distance between two points D = 1km /cos(theta)
|
---|
568 | // Density along y within this km: f = (x[i+1]-x[i])/1km * dy
|
---|
569 | // Integral of this density f = (x[i+1]-x[i])/1km * (y[i+1]-y[i])
|
---|
570 | // f(h) = int [ (c1-c0)/1km*(h-h0)*dh + c0 ] dh
|
---|
571 | // = (c-co)*(h-h0)
|
---|
572 |
|
---|
573 | Double_t rc = 0;
|
---|
574 | int i;
|
---|
575 | for (i=0; i<49; i++)
|
---|
576 | if (i>=2 && i+1<height/1e5) // cm -> km
|
---|
577 | rc += oz_conc[i] * 1e5/cos(theta);
|
---|
578 |
|
---|
579 | rc -= oz_conc[2]*0.2*1e5/cos(theta);
|
---|
580 | rc += oz_conc[i+1]*fmod(height/1e5,1)*1e5/cos(theta);
|
---|
581 |
|
---|
582 | return rc;
|
---|
583 | }
|
---|
584 | */
|
---|
585 |
|
---|
586 | Double_t CalcOzoneAbsorption(Double_t h, Double_t wavelength, Double_t theta) const
|
---|
587 | {
|
---|
588 | if (!fAbsCoeffOzone)
|
---|
589 | return 1;
|
---|
590 |
|
---|
591 | //******* Ozone absorption *******
|
---|
592 | if (h > 50.e5)
|
---|
593 | h = 50.e5;
|
---|
594 |
|
---|
595 | // Vigroux Ozone absorption coefficient a.s.l. through interpolation:
|
---|
596 | //const float oz_vigroux[15]= {1.06e2, 1.01e1, 8.98e-1, 6.40e-2, 1.80e-3, 0, 0, 3.50e-3, 3.45e-2, 9.20e-2, 1.32e-1, 6.20e-2, 2.30e-2, 1.00e-2, 0.00};
|
---|
597 | //const Double_t beta0 = getbeta(wavelength, oz_vigroux);
|
---|
598 | const Double_t beta0 = GetBeta(wavelength, *fAbsCoeffOzone);
|
---|
599 |
|
---|
600 | // Now use the pre-calculated values of the path integral
|
---|
601 | // for h and theta
|
---|
602 | const UInt_t H = TMath::Nint(h/1e4);
|
---|
603 | const UInt_t T = TMath::Min(89, TMath::Nint(theta/STEPTHETA));
|
---|
604 |
|
---|
605 | const Double_t path = ozone_path[H][T];
|
---|
606 |
|
---|
607 | return TMath::Exp(-beta0*path);
|
---|
608 | }
|
---|
609 |
|
---|
610 | Double_t CalcAerosolAbsorption(Double_t h, Double_t wavelength, Double_t theta) const
|
---|
611 | {
|
---|
612 | if (!fAbsCoeffAerosols)
|
---|
613 | return 1;
|
---|
614 |
|
---|
615 | //******* Mie (aerosol) *******
|
---|
616 | if (h > 30.e5)
|
---|
617 | h = 30.e5;
|
---|
618 |
|
---|
619 | // Mie absorption coefficient at sea level through interpolation
|
---|
620 | //const float aero_betap[15] = {0.27, 0.26, 0.25, 0.24, 0.24, 0.23, 0.20, 0.180, 0.167, 0.158, 0.150, 0.142, 0.135, 0.127, 0.120};
|
---|
621 | //const Double_t beta0 = getbeta(wavelength, aero_betap);
|
---|
622 | const Double_t beta0 = GetBeta(wavelength, *fAbsCoeffAerosols);
|
---|
623 |
|
---|
624 | // Now use the pre-calculated values of the path integral
|
---|
625 | // for h and theta
|
---|
626 | const UInt_t H = TMath::Nint(h/1e4);
|
---|
627 | const UInt_t T = TMath::Min(89, TMath::Nint(theta/STEPTHETA));
|
---|
628 |
|
---|
629 | const Double_t path = aerosol_path[H][T];
|
---|
630 |
|
---|
631 | return TMath::Exp(-beta0*path);
|
---|
632 | }
|
---|
633 |
|
---|
634 | Double_t GetTransmission(const MPhotonData &ph) const
|
---|
635 | {
|
---|
636 | const Double_t wavelength = ph.GetWavelength();
|
---|
637 | const Double_t height = ph.GetProductionHeight();
|
---|
638 |
|
---|
639 | // Reduce the necessary number of floating point operations
|
---|
640 | // by storing the intermediate results
|
---|
641 | const Double_t sin2 = ph.GetSinW2();
|
---|
642 | const Double_t cost = TMath::Sqrt(1-sin2);
|
---|
643 | const Double_t theta = TMath::ACos(cost);
|
---|
644 |
|
---|
645 | // Path from production height to obslevel
|
---|
646 | const Double_t z = height-fObsLevel;
|
---|
647 |
|
---|
648 | // Distance of emission point to incident point on ground
|
---|
649 | const Double_t d = z/cost;
|
---|
650 |
|
---|
651 | // Avoid problems if photon is very close to telescope:
|
---|
652 | if (TMath::Abs(d)<1)
|
---|
653 | return 1;
|
---|
654 |
|
---|
655 | // Earth radius plus observation height (distance of telescope
|
---|
656 | // from earth center)
|
---|
657 | const Double_t H = R() + fObsLevel;
|
---|
658 |
|
---|
659 | // We calculate h, the true height a.s.l.
|
---|
660 | // of the photon emission point in cm
|
---|
661 | const Double_t h = TMath::Sqrt(H*H + d*d + 2*H*z) - R();
|
---|
662 |
|
---|
663 | //**** Rayleigh scattering: *****
|
---|
664 | const Double_t T_Ray = CalcTransmission(h, wavelength, sin2);
|
---|
665 | if (T_Ray<0)
|
---|
666 | return 0;
|
---|
667 |
|
---|
668 | //****** Ozone absorption: ******
|
---|
669 | const Double_t T_Oz = CalcOzoneAbsorption(h, wavelength, theta);
|
---|
670 |
|
---|
671 | //******** Mie (aerosol) ********
|
---|
672 | const Double_t T_Mie = CalcAerosolAbsorption(h, wavelength, theta);
|
---|
673 |
|
---|
674 | // FIXME: What if I wanna display these values?
|
---|
675 |
|
---|
676 | // Calculate final transmission coefficient
|
---|
677 | return T_Ray * T_Oz * T_Mie;
|
---|
678 | }
|
---|
679 | };
|
---|
680 |
|
---|
681 | const Double_t MAtmosphere::STEPTHETA = 1.74533e-2; // aprox. 1 degree
|
---|
682 |
|
---|
683 | const Double_t MAtmRayleigh::fgMeanFreePath = 2970;
|
---|
684 |
|
---|
685 | const Double_t MAtmosphere::aero_n[31] = {200, 87, 38, 16, 7.2, 3.1, 1.1, 0.4, 0.14, 5.0e-2, 2.6e-2, 2.3e-2, 2.1e-2, 2.3e-2, 2.5e-2, 4.1e-2, 6.7e-2, 7.3e-2, 8.0e-2, 9.0e-2, 8.6e-2, 8.2e-2, 8.0e-2, 7.6e-2, 5.2e-2, 3.6e-2, 2.5e-2, 2.4e-2, 2.2e-2, 2.0e-2, 1.9e-2};
|
---|
686 |
|
---|
687 | const Double_t MAtmosphere::oz_conc[51]={0.3556603E-02, 0.3264150E-02, 0.2933961E-02, 0.2499999E-02, 0.2264150E-02, 0.2207546E-02, 0.2160377E-02, 0.2226414E-02, 0.2283018E-02, 0.2811320E-02, 0.3499999E-02, 0.4603772E-02, 0.6207545E-02, 0.8452828E-02, 0.9528299E-02, 0.9905657E-02, 0.1028302E-01, 0.1113207E-01, 0.1216981E-01, 0.1424528E-01, 0.1641509E-01, 0.1839622E-01, 0.1971697E-01, 0.1981131E-01, 0.1933962E-01, 0.1801886E-01, 0.1632075E-01, 0.1405660E-01, 0.1226415E-01, 0.1066037E-01, 0.9028300E-02, 0.7933960E-02, 0.6830187E-02, 0.5820753E-02, 0.4830188E-02, 0.4311319E-02, 0.3613206E-02, 0.3018867E-02, 0.2528301E-02, 0.2169811E-02, 0.1858490E-02, 0.1518867E-02, 0.1188679E-02, 0.9301884E-03, 0.7443394E-03, 0.5764149E-03, 0.4462263E-03, 0.3528301E-03, 0.2792452E-03, 0.2226415E-03, 0.1858490E-03};
|
---|
688 |
|
---|
689 | // ==========================================================================
|
---|
690 |
|
---|
691 | // --------------------------------------------------------------------------
|
---|
692 | //
|
---|
693 | // Default Constructor.
|
---|
694 | //
|
---|
695 | MSimAtmosphere::MSimAtmosphere(const char* name, const char *title)
|
---|
696 | : fEvt(0), fAtmosphere(0)
|
---|
697 | {
|
---|
698 | fName = name ? name : "MSimAtmosphere";
|
---|
699 | fTitle = title ? title : "Simulate the wavelength and height-dependant atmpsheric absorption";
|
---|
700 |
|
---|
701 | fAtmosphere = new MAtmosphere;
|
---|
702 | }
|
---|
703 |
|
---|
704 | // --------------------------------------------------------------------------
|
---|
705 | //
|
---|
706 | // Calls Clear()
|
---|
707 | //
|
---|
708 | MSimAtmosphere::~MSimAtmosphere()
|
---|
709 | {
|
---|
710 | delete fAtmosphere;
|
---|
711 | }
|
---|
712 |
|
---|
713 | // --------------------------------------------------------------------------
|
---|
714 | //
|
---|
715 | // Search for the needed parameter containers. Read spline from file
|
---|
716 | // calling ReadFile();
|
---|
717 | //
|
---|
718 | Int_t MSimAtmosphere::PreProcess(MParList *pList)
|
---|
719 | {
|
---|
720 | fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
|
---|
721 | if (!fEvt)
|
---|
722 | {
|
---|
723 | *fLog << err << "MPhotonEvent not found... aborting." << endl;
|
---|
724 | return kFALSE;
|
---|
725 | }
|
---|
726 |
|
---|
727 |
|
---|
728 | return kTRUE;
|
---|
729 | }
|
---|
730 |
|
---|
731 | // --------------------------------------------------------------------------
|
---|
732 | //
|
---|
733 | Bool_t MSimAtmosphere::ReInit(MParList *pList)
|
---|
734 | {
|
---|
735 | MCorsikaRunHeader *h = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
|
---|
736 | if (!h)
|
---|
737 | {
|
---|
738 | *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
|
---|
739 | return kFALSE;
|
---|
740 | }
|
---|
741 |
|
---|
742 | //if (fRunHeader->Has(MCorsikaRunHeader::kRefraction))
|
---|
743 | // *fLog << inf << "Atmospheric refraction already applied in Corsika... skipping our own." << endl;
|
---|
744 |
|
---|
745 | // FIXME: Check wavelength range
|
---|
746 |
|
---|
747 | /*
|
---|
748 | if (h->GetWavelengthMin()<fSpline->GetXmin())
|
---|
749 | *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds lower bound of spline." << endl;
|
---|
750 |
|
---|
751 | if (h->GetWavelengthMax()>fSpline->GetXmax())
|
---|
752 | *fLog << warn << "WARNING - Upper bound of wavelength bandwidth exceeds upper bound of spline." << endl;
|
---|
753 | */
|
---|
754 |
|
---|
755 | fAtmosphere->Init(*h);
|
---|
756 |
|
---|
757 | if (!fAtmosphere->IsAllValid())
|
---|
758 | {
|
---|
759 | *fLog << err << "ERROR - Something with the atmoshere's initialization went wrong!" << endl;
|
---|
760 | return kFALSE;
|
---|
761 | }
|
---|
762 |
|
---|
763 | if (h->GetWavelengthMin()<fAtmosphere->GetWavelengthMin())
|
---|
764 | *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds valid range of atmosphere." << endl;
|
---|
765 |
|
---|
766 | if (h->GetWavelengthMax()>fAtmosphere->GetWavelengthMax())
|
---|
767 | *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds valid range of atmosphere." << endl;
|
---|
768 |
|
---|
769 | if (!h->Has(MCorsikaRunHeader::kAtmext))
|
---|
770 | *fLog << warn << "WARNING - ATMEXT option not used for Corsika data." << endl;
|
---|
771 |
|
---|
772 | if (!h->Has(MCorsikaRunHeader::kRefraction))
|
---|
773 | *fLog << warn << "WARNING - Refraction calculation disabled for Corsika data." << endl;
|
---|
774 |
|
---|
775 | return kTRUE;
|
---|
776 | }
|
---|
777 |
|
---|
778 | // --------------------------------------------------------------------------
|
---|
779 | //
|
---|
780 | Int_t MSimAtmosphere::Process()
|
---|
781 | {
|
---|
782 | // Get the number of photons in the list
|
---|
783 | const Int_t num = fEvt->GetNumPhotons();
|
---|
784 |
|
---|
785 | // FIMXE: Add checks for
|
---|
786 | // * upgoing particles
|
---|
787 | // * Can we take the full length until the camera into account?
|
---|
788 |
|
---|
789 | // Counter for number of total and final events
|
---|
790 | Int_t cnt = 0;
|
---|
791 | for (Int_t i=0; i<num; i++)
|
---|
792 | {
|
---|
793 | // Get i-th photon from the list
|
---|
794 | const MPhotonData &ph = (*fEvt)[i];
|
---|
795 |
|
---|
796 | // Get atmospheric transmission for this photon
|
---|
797 | const Double_t eff = fAtmosphere->GetTransmission(ph);
|
---|
798 |
|
---|
799 | // Get a random value between 0 and 1 to determine whether the photon will survive
|
---|
800 | // gRandom->Rndm() = [0;1[
|
---|
801 | if (gRandom->Rndm()>=eff)
|
---|
802 | continue;
|
---|
803 |
|
---|
804 | // Copy the surviving events bakc in the list
|
---|
805 | (*fEvt)[cnt++] = ph;
|
---|
806 | }
|
---|
807 |
|
---|
808 | // Now we shrink the array to the number of new entries.
|
---|
809 | fEvt->Shrink(cnt);
|
---|
810 |
|
---|
811 | return kTRUE;
|
---|
812 | }
|
---|
813 |
|
---|
814 | /*
|
---|
815 | Int_t MSimWavelength::Process()
|
---|
816 | {
|
---|
817 | // Get the number of photons in the list
|
---|
818 | const Int_t num = fEvt->GetNumPhotons();
|
---|
819 |
|
---|
820 | // FIMXE: Add checks for
|
---|
821 | // * upgoing particles
|
---|
822 | // * wavelength range
|
---|
823 | // * check if corsika atmosphere is switched on
|
---|
824 | // * Can we take the full length until the camera into account?
|
---|
825 |
|
---|
826 | // Counter for number of total and final events
|
---|
827 | Int_t cnt = 0;
|
---|
828 | for (Int_t i=0; i<num; i++)
|
---|
829 | {
|
---|
830 | // Get i-th photon from the list
|
---|
831 | MPhotonData &ph = (*fEvt)[i];
|
---|
832 |
|
---|
833 | const Double_t min = fRunHeader->GetWavelengthMin(); // WAVLGL
|
---|
834 | const Double_t max = fRunHeader->GetWavelengthMax(); // WAVLGU
|
---|
835 | const Double_t f = (max-min)/max;
|
---|
836 |
|
---|
837 | // WAVELENGTH = 1. / (1/min - RD(1)/(min*max/(max-min)))
|
---|
838 |
|
---|
839 |
|
---|
840 | ph.SetWavelength(TMath::Nint(min / (1. - gRandom->Rndm()*f)));
|
---|
841 | }
|
---|
842 |
|
---|
843 | return kTRUE;
|
---|
844 | }
|
---|
845 | */
|
---|
846 |
|
---|
847 | // --------------------------------------------------------------------------
|
---|
848 | //
|
---|
849 | // FileName: reflectivity.txt
|
---|
850 | // UseTheta: No
|
---|
851 | //
|
---|
852 | Int_t MSimAtmosphere::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
|
---|
853 | {
|
---|
854 | Bool_t rc = kFALSE;
|
---|
855 | /*
|
---|
856 | if (IsEnvDefined(env, prefix, "FileName", print))
|
---|
857 | {
|
---|
858 | rc = kTRUE;
|
---|
859 | SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
|
---|
860 | }
|
---|
861 |
|
---|
862 | if (IsEnvDefined(env, prefix, "UseTheta", print))
|
---|
863 | {
|
---|
864 | rc = kTRUE;
|
---|
865 | SetUseTheta(GetEnvValue(env, prefix, "UseTheta", fUseTheta));
|
---|
866 | }
|
---|
867 | */
|
---|
868 | return rc;
|
---|
869 | }
|
---|