Index: trunk/Mars/msim/MAtmosphere.cc =================================================================== --- trunk/Mars/msim/MAtmosphere.cc (revision 19763) +++ trunk/Mars/msim/MAtmosphere.cc (revision 19763) @@ -0,0 +1,628 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of CheObs, the Modular Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appears in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 1/2009 +! +! Copyright: CheObs Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSimAtmosphere +// +// Task to calculate wavelength or incident angle dependent absorption +// +// Input Containers: +// MPhotonEvent +// MCorsikaRunHeader +// +// Output Containers: +// MPhotonEvent +// +////////////////////////////////////////////////////////////////////////////// +#include "MAtmosphere.h" + +#include + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCorsikaRunHeader.h" +#include "MPhotonEvent.h" +#include "MPhotonData.h" + +ClassImp(MAtmosphere); +ClassImp(MAtmRayleigh); + +using namespace std; + +// ========================================================================== +// +// January 2002, A. Moralejo: We now precalculate the slant paths for the +// aerosol and Ozone vertical profiles, and then do an interpolation in +// wavelength for every photon to get the optical depths. The parameters +// used, defined below, have been taken from "Atmospheric Optics", by +// L. Elterman and R.B. Toolin, chapter 7 of the "Handbook of geophysics +// and Space environments". (S.L. Valley, editor). McGraw-Hill, NY 1965. +// +// WARNING: the Mie scattering and the Ozone absorption are implemented +// to work only with photons produced at a height a.s.l larger than the +// observation level. So this is not expected to work well for simulating +// the telescope pointing at theta > 90 deg (for instance for neutrino +// studies. Rayleigh scattering works even for light coming from below. +// +// Fixed bugs (of small influence) in Mie absorption implementation: there +// were errors in the optical depths table, as well as a confusion: +// height a.s.l. was used as if it was height above the telescope level. +// The latter error was also present in the Ozone aborption implementation. +// +// On the other hand, now we have the tables AE_ABI and OZ_ABI with optical +// depths between sea level and a height h (before it was between 2km a.s.l +// and a height h). So that we can simulate also in the future a different +// observation level. +// +// AM: WARNING: IF VERY LARGE zenith angle simulations are to be done (say +// above 85 degrees, for neutrino primaries or any other purpose) this code +// will have to be adapted accordingly and checked, since in principle it has +// been written and tested to simulate the absorption of Cherenkov photons +// arriving at the telescope from above. +// +// AM: WARNING 2: not to be used for wavelengths outside the range 250-700 nm +// +// January 2003, Abelardo Moralejo: found error in Ozone absorption treatment. +// At large zenith angles, the air mass used was the one calculated for +// Rayleigh scattering, but since the Ozone distribution is rather different +// from the global distribution of air molecules, this is not a good +// approximation. Now I have left in this code only the Rayleigh scattering, +// and moved to atm.c the Mie scattering and Ozone absorption calculated in +// a better way. +// +// A. Moralejo, January 2003: added some parameters for Mie scattering +// and Ozone absorption derived from the clear standard atmosphere model +// in "Atmospheric Optics", by L. Elterman and R.B. Toolin, chapter 7 of +// the "Handbook of geophysics and Space environments". S.L. Valley, +// editor. McGraw-Hill, NY 1965. +// +// AM, Jan 2003: Changed the meaning of the argument height: now it is the +// true height above sea level at which a photon has been emitted, before +// it was the height given by Corsika, measured in the vertical of the +// observer and not in the vertical of the emitting particle. +// +// +// MAGIC-Winter and MAGIC-Summer by M. Haffke, +// parametrizing profiles obtained with MSIS: +// http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm +// +// +// The MAGIC-Winter and MAGIC-Summer parametrisations reproduce the MSIS +// profiles for the 3 atmospheric layers from 0 up to 40 km height. Beyond +// that height, it was not possible to achieve a good fit, but the amount +// of residual atmosphere is so small that light absorption would be +// negligible anyway. Showers develop well below 40 km. +// +// +// The mass overburden is given by T = AATM + BATM * exp(-h/CATM) +// The last layer of the US standard atmosphere (in which T varies +// linearly with h) is above 100 km and has not been included here +// because it should not matter. +// + +const Double_t MAtmosphere::STEPTHETA = 1.74533e-2; // aprox. 1 degree + +const Double_t MAtmRayleigh::fgMeanFreePath = 2970; + +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}; + +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}; + +MAtmRayleigh::MAtmRayleigh() : fR(MCorsikaRunHeader::fgEarthRadius), + fHeight{0, 300000, 1e+06, 5e+06, 1e+07}, + //fAtmA{-144.838, -124.071, 0.360027, -0.000824761, 0.00221589}, + fAtmB{1192.34, 1173.98, 1412.08, 810.682, 1}, + fAtmC{994186, 967530, 636143, 735640, 4.92961e9}, + fObsLevel(-1) +{ +} + +// -------------------------------------------------------------------------- +// +// Precalcalculate the integrals from the observer level to the next +// atmpsheric layer below including the lower boundary. Thus a +// correct calculation is reduced to the calculation of the upper +// boundary. +// +// fRho[0] = B0; +// fRho[1] = B0-A0 + B1; +// fRho[2] = B0-A0 + B1-A1 + B2; +// fRho[3] = B0-A0 + B1-A1 + B2+A2 + B3; +// fRho[4] = B0-A0 + B1-A1 + B2+A2 + B3 - A3; +// +void MAtmRayleigh::PreCalcRho() +{ + // Limits (height in cm) of the four layers in which + // atmosphere is parametrized. + // This is a stupid trick giving 0 for the integrals below + // the observer + + // FIXME: Could be replaced by 0, AtmLay[0]-fAtmLay[3] + + const Double_t h[5] = + { + fObsLevel, // fObsLevel, // 0km + TMath::Max(fObsLevel, 7.75e5), // TMath::Max(fObsLevel, 4e5), // 4km + 16.5e5, // 10e5, // 10km + 50.0e5, // 40e5, // 40km + 105.0e5, // 100e5 // 100km + }; + + memcpy(fHeight, h, sizeof(Double_t)*5); + + fRho[0] = 0; + for (int i=0; i<4; i++) + { + const Double_t b = fAtmB[i]; + const Double_t c = fAtmC[i]; + + const Double_t h1 = h[i+1]; + const Double_t h0 = h[i]; + + const Double_t B = b*TMath::Exp(-h0/c); + const Double_t A = b*TMath::Exp(-h1/c); + + // Calculate rho for the i-th layer from the lower + // to the higher layer boundary. + // If height is within the layer only calculate up to height. + fRho[i] += B; + fRho[i+1] = fRho[i] - A; + } +} + +void MAtmRayleigh::Init(Double_t obs, const Float_t *atmb, const Float_t *atmc) +{ + // Observation level above earth radius + fObsLevel = obs; + + //memcpy(fAtmA, (Float_t*)h.GetAtmosphericCoeffA(), sizeof(Float_t)*4); + if (atmb) + memcpy(fAtmB, atmb, sizeof(Float_t)*5); + if (atmc) + memcpy(fAtmC, atmc, sizeof(Float_t)*5); + + PreCalcRho(); +} + +// Init an atmosphere from the data stored in MCorsikaRunHeader +// This initialized fObsLevel, fR, fAtmB and fAtmC and +// PreCalcRho +void MAtmRayleigh::Init(const MCorsikaRunHeader &h) +{ + // Use earth radius as defined in Corsika + fR = h.EarthRadius(); + + Init(h.GetObsLevel(), h.GetAtmosphericCoeffB(), h.GetAtmosphericCoeffC()); +} + +// Return the vertical thickness between the observer and height. +// Therefor the integral of the layers below (including the lower +// boudary) Have been precalculated by PreCalcRho +Double_t MAtmRayleigh::GetVerticalThickness(Double_t height) const +{ + // FIXME: We could store the start point above obs-level + // (Does this really gain anything?) + Int_t i=0; + while (i<4 && height>fHeight[i+1]) + i++; + + const Double_t b = fAtmB[i]; + const Double_t c = fAtmC[i]; + + return fRho[i] - b*TMath::Exp(-height/c); +} + +/* +// The "orginal" code for the vertical thickness +Double_t GetVerticalThickness(Double_t obslev, Double_t height) const +{ + // This is a C++-version of the original code from attenu.c + + // Limits (height in cm) of the four layers in which atmosphere is parametrized: + const double lahg[5] = + { + obslev, + TMath::Max(obslev, 4e5), + 1.0e6, + 4.0e6, + 1.0e7 + }; + + Double_t Rho_Tot = 0.0; + for (int i=0; i<4; i++) + { + const Double_t b = fAtmB[i]; + const Double_t c = fAtmC[i]; + + const Double_t h0 = TMath::Min(height, lahg[i+1]); + const Double_t h1 = lahg[i]; + + // Calculate rho for the i-th layer from the lower + // to the higher layer boundary. + // If height is within the layer only calculate up to height. + Rho_Tot += b*(exp(-h1/c) - exp(-h0/c)); + + if (lahg[i+1] > height) + break; + } + + return Rho_Tot; +} +*/ + +Double_t MAtmRayleigh::CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const +{ + // sin2: sin(theta)^2 + // height is the true height a.s.l. + + // LARGE ZENITH ANGLE FACTOR (AIR MASS FACTOR): + // Air mass factor "airmass" calculated using a one-exponential + // density profile for the atmosphere, + // + // rho = rho_0 exp(-height/hscale) with hscale = 7.4 km + // + // The air mass factor is defined as I(theta)/I(0), the ratio of + // the optical paths I (in g/cm2) traversed between two given + // heights, at theta and at 0 deg z.a. + + const Double_t H = height-fObsLevel; + const Double_t h = 2*H; + + // Scale-height (cm) for Exponential density profile + const Double_t hscale = 7.4e5; + const Double_t f = 2*hscale; + + // Precalc R*cos(theta)^2 (FIXME: Is ph.GetCosW2 more precise?) + const Double_t Rcos2 = fR * (1-sin2); // cos2 = 1 - sin2 + + const Double_t x1 = TMath::Sqrt((Rcos2 ) / f); + const Double_t x2 = TMath::Sqrt((Rcos2 + 2*h) / f); + const Double_t x3 = TMath::Sqrt((fR ) / f); + const Double_t x4 = TMath::Sqrt((fR + 2*h) / f); + + // Return a -1 transmittance in the case the photon comes + // exactly from the observation level, because in that case the + // "air mass factor" would become infinity and the calculation + // is not valid! + if (fabs(x3-x4) < 1.e-10) + return -1.; + + const Double_t e12 = erfc(x1) - erfc(x2); + const Double_t e34 = erfc(x3) - erfc(x4); + + const Double_t airmass = TMath::Exp(-fR*sin2 / f) * e12/e34; + + // Calculate the traversed "vertical thickness" of air using the + // US Standard atmosphere: + const Double_t Rho_tot = GetVerticalThickness(/*fObsLevel,*/ height); + + // We now convert from "vertical thickness" to "slanted thickness" + // traversed by the photon on its way to the telescope, simply + // multiplying by the air mass factor m: + const Double_t Rho_Fi = airmass * Rho_tot; + + // Finally we calculate the transmission coefficient for the Rayleigh + // scattering: + // AM Dec 2002, introduced ABS below to account (in the future) for + // possible photons coming from below the observation level. + + const Double_t wl = 400./wavelength; + + // Mean free path for scattering Rayleigh XR (g/cm^2) + return TMath::Exp(-TMath::Abs(Rho_Fi/fgMeanFreePath)*wl*wl*wl*wl); +} + +// ========================================================================== + +// Interpolate the graph at wavelength +Double_t MAtmosphere::GetBeta(Double_t wavelength, const TGraph &g) const +{ + // FIXME: This might not be the fastest because range + // checks are done for each call! + return g.GetN()==0 ? 0 : g.Eval(wavelength)*1e-5; // from km^-1 to cm^-1 +/* + // Linear interpolation + // (FIXME: Is it faster to be replaced with a binary search?) + // ( This might be faster because we have more photons + // with smaller wavelengths) + //int index; + //for (index = 1; index GetX()[0], fAbsCoeffAerosols->GetX()[0]) : -1; +} + +Float_t MAtmosphere::GetWavelengthMax() const +{ + return fAbsCoeffOzone && fAbsCoeffAerosols ? TMath::Min(fAbsCoeffOzone->GetX()[fAbsCoeffOzone->GetN()-1], fAbsCoeffAerosols->GetX()[fAbsCoeffAerosols->GetN()-1]) : -1; +} + +Bool_t MAtmosphere::HasValidOzone() const +{ + return fAbsCoeffOzone && fAbsCoeffOzone->GetN()>0; +} + +Bool_t MAtmosphere::HasValidAerosol() const +{ + return fAbsCoeffAerosols && fAbsCoeffAerosols->GetN()>0; +} + +void MAtmosphere::PreCalcOzone() +{ + // It follows a precalculation of the slant path integrals we need + // for the estimate of the Mie scattering and Ozone absorption: + Double_t dh = 1.e3; + //const Double_t STEPTHETA = 1.74533e-2; // aprox. 1 degree + + // Ozone absorption + for (Int_t j = 0; j < 90; j++) + { + const Double_t theta = j * STEPTHETA; + const Double_t sin2 = sin(theta)*sin(theta); + const Double_t H = R()+fObsLevel; + + Double_t path_slant = 0; + for (Double_t h = fObsLevel; h <= 50e5; h += dh) + { + // h is the true height vertical above ground + if (fmod(h,1e4) == 0) + ozone_path[(int)(h/1e4)][j] = path_slant; + + const Double_t km = h/1e5; + const Int_t i = TMath::FloorNint(km); + const Double_t l = R()+h; + + const Double_t L = TMath::Sqrt(l*l - H*H * sin2); + const Double_t f = dh * l / L; + + // Linear interpolation at h/1e5 + Double_t interpol = oz_conc[i] + fmod(km, 1) * (oz_conc[i+1]-oz_conc[i]); + + path_slant += f * interpol; + } + } +} + +void MAtmosphere::PreCalcAerosol() +{ + // It follows a precalculation of the slant path integrals we need + // for the estimate of the Mie scattering and Ozone absorption: + Double_t dh = 1.e3; + //const Double_t STEPTHETA = 1.74533e-2; // aprox. 1 degree + + /* Mie (aerosol): */ + for (Int_t j = 0; j < 90; j++) + { + const Double_t theta = j * STEPTHETA; + const Double_t sin2 = sin(theta)*sin(theta); + const Double_t H = R()+fObsLevel; + + Double_t path_slant = 0; + for (Double_t h = fObsLevel; h <= 30e5; h += dh) + { + // h is the true height vertical above ground + if (fmod(h,1e4) == 0) + aerosol_path[(int)(h/1e4)][j] = path_slant; + + const Double_t km = h/1e5; + const Int_t i = TMath::FloorNint(km); + const Double_t l = R()+h; + + const Double_t L = TMath::Sqrt(l*l - H*H * sin2); + const Double_t f = dh * l / L; + + // Linear interpolation at h/1e5 + Double_t interpol = aero_n[i] + fmod(km, 1)*(aero_n[i+1]-aero_n[i]); + + path_slant += f * interpol/aero_n[0]; // aero_n [km^-1] + } + } +} + +Bool_t MAtmosphere::InitOzone(const TString name) +{ + if (!name.IsNull()) + { + if (fAbsCoeffOzone) + delete fAbsCoeffOzone; + + fAbsCoeffOzone = new TGraph(name); + fAbsCoeffOzone->Sort(); + } + + if (!HasValidAerosol()) + return kFALSE; + + if (IsValid()) + PreCalcOzone(); + + return kTRUE; +} + +Bool_t MAtmosphere::InitAerosols(const TString name) +{ + if (!name.IsNull()) + { + if (fAbsCoeffAerosols) + delete fAbsCoeffAerosols; + + fAbsCoeffAerosols = new TGraph(name); + fAbsCoeffAerosols->Sort(); + } + + if (!HasValidAerosol()) + return kFALSE; + + if (IsValid()) + PreCalcAerosol(); + + return kTRUE; +} + +/* +Double_t GetOz(Double_t height, Double_t theta) const +{ + // Distance between two points D = 1km /cos(theta) + // Density along y within this km: f = (x[i+1]-x[i])/1km * dy + // Integral of this density f = (x[i+1]-x[i])/1km * (y[i+1]-y[i]) + // f(h) = int [ (c1-c0)/1km*(h-h0)*dh + c0 ] dh + // = (c-co)*(h-h0) + + Double_t rc = 0; + int i; + for (i=0; i<49; i++) + if (i>=2 && i+1 km + rc += oz_conc[i] * 1e5/cos(theta); + + rc -= oz_conc[2]*0.2*1e5/cos(theta); + rc += oz_conc[i+1]*fmod(height/1e5,1)*1e5/cos(theta); + + return rc; +} +*/ + +Double_t MAtmosphere::CalcOzoneAbsorption(Double_t h, Double_t wavelength, Double_t theta) const +{ + if (!fAbsCoeffOzone) + return 1; + + //******* Ozone absorption ******* + if (h > 50.e5) + h = 50.e5; + + // Vigroux Ozone absorption coefficient a.s.l. through interpolation: + //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}; + //const Double_t beta0 = getbeta(wavelength, oz_vigroux); + const Double_t beta0 = GetBeta(wavelength, *fAbsCoeffOzone); + + // Now use the pre-calculated values of the path integral + // for h and theta + const UInt_t H = TMath::Nint(h/1e4); + const UInt_t T = TMath::Min(89, TMath::Nint(theta/STEPTHETA)); + + const Double_t path = ozone_path[H][T]; + + return TMath::Exp(-beta0*path); +} + +Double_t MAtmosphere::CalcAerosolAbsorption(Double_t h, Double_t wavelength, Double_t theta) const +{ + if (!fAbsCoeffAerosols) + return 1; + + //******* Mie (aerosol) ******* + if (h > 30.e5) + h = 30.e5; + + //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}; + //const Double_t beta0 = getbeta(wavelength, aero_betap); + const Double_t beta0 = GetBeta(wavelength, *fAbsCoeffAerosols); + + // Now use the pre-calculated values of the path integral + // for h and theta + const UInt_t H = TMath::Nint(h/1e4); + const UInt_t T = TMath::Min(89, TMath::Nint(theta/STEPTHETA)); + + + const Double_t path = aerosol_path[H][T]; + + return TMath::Exp(-beta0*path); +} + +Double_t MAtmosphere::GetTransmission(const MPhotonData &ph) const +{ + const Double_t wavelength = ph.GetWavelength(); + const Double_t height = ph.GetProductionHeight(); + + // Reduce the necessary number of floating point operations + // by storing the intermediate results + const Double_t sin2 = ph.GetSinW2(); + const Double_t cost = TMath::Sqrt(1-sin2); + const Double_t theta = TMath::ACos(cost); + + // Path from production height to obslevel + const Double_t z = height-fObsLevel; + + // Distance of emission point to incident point on ground + const Double_t d = z/cost; + + // Avoid problems if photon is very close to telescope: + if (TMath::Abs(d)<1) + return 1; + + // Earth radius plus observation height (distance of telescope + // from earth center) + const Double_t H = R() + fObsLevel; + + // We calculate h, the true height a.s.l. + // of the photon emission point in cm + const Double_t h = TMath::Sqrt(H*H + d*d + 2*H*z) - R(); + + //**** Rayleigh scattering: ***** + const Double_t T_Ray = CalcTransmission(h, wavelength, sin2); + if (T_Ray<0) + return 0; + + //****** Ozone absorption: ****** + const Double_t T_Oz = CalcOzoneAbsorption(h, wavelength, theta); + + //******** Mie (aerosol) ******** + const Double_t T_Mie = CalcAerosolAbsorption(h, wavelength, theta); + + // FIXME: What if I wanna display these values? + + // Calculate final transmission coefficient + return T_Ray * T_Oz * T_Mie; +} + Index: trunk/Mars/msim/MAtmosphere.h =================================================================== --- trunk/Mars/msim/MAtmosphere.h (revision 19763) +++ trunk/Mars/msim/MAtmosphere.h (revision 19763) @@ -0,0 +1,137 @@ +#ifndef MARS_MAtmosphere +#define MARS_MAtmosphere + +#include + +class TGraph; +class MPhotonData; +class MCorsikaRunHeader; + +class MAtmRayleigh +{ +private: + static const Double_t fgMeanFreePath; // [g/cm^2] Mean free path for scattering Rayleigh XR + static const Double_t fgEarthRadiush; // [cm] Default Earth Radius + + Double_t fR; // [cm] Earth radius to be used + + Double_t fHeight[5]; // Layer boundaries (atmospheric layer) + + // Parameters of the different atmospheres. We use the same parametrization + // shape as in Corsika atmospheric models (see Corsika manual, appendix D). + // The values here can be/are obtained from the Corsika output + //Float_t fAtmA[5]; // The index refers to the atmospheric layer (starting from sea level and going upwards) + Float_t fAtmB[5]; // The index refers to the atmospheric layer (starting from sea level and going upwards) + Float_t fAtmC[5]; // The index refers to the atmospheric layer (starting from sea level and going upwards) + + Double_t fRho[5]; // Precalculated integrals for rayleigh scatterning + + void PreCalcRho(); + +protected: + Double_t fObsLevel; // [cm] observation level a.s.l. + +public: + // Default constructor + MAtmRayleigh(); + + // Init an atmosphere from the data stored in MCorsikaRunHeader + MAtmRayleigh(const MCorsikaRunHeader &h) + { + Init(h); + } + + // Check if the ovservation level has been correctly initialized + // Used as a marker for correct initialization + Bool_t IsValid() const { return fObsLevel>=0; } + + // Get the Earth radius to be used + Double_t R() const { return fR; } + + void Init(Double_t obs, const Float_t *atmb=0, const Float_t *atmc=0); + void Init(const MCorsikaRunHeader &h); + + Double_t GetVerticalThickness(Double_t height) const; + Double_t CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const; +}; + +// ========================================================================== + +class MAtmosphere : public MAtmRayleigh +{ +private: + static const Double_t STEPTHETA; // aprox. 1 degree + + // Aerosol number density for 31 heights a.s.l., from 0 to 30 km, + // in 1 km steps (units: cm^-3) + static const Double_t aero_n[31]; + + // Ozone concentration for 51 heights a.s.l., from 0 to 50 km, + // in 1 km steps (units: cm/km) + static const Double_t oz_conc[51]; + + // aerosol_path contains the path integrals for the aerosol number + // density (relative to the number density at sea level) between the + // observation level and a height h for different zenith angles. The + // first index indicate height above sea level in units of 100m, the + // second is the zenith angle in degrees. + float aerosol_path[301][90]; + + // ozone_path contains the path integrals for the ozone concentration + // between the observation level and a height h for different zenith + // angles. The first index indicate height above sea level in units + // of 100m, the second is the zenith angle in degrees. + float ozone_path[501][90]; + + // Interpolate the graph at wavelength + Double_t GetBeta(Double_t wavelength, const TGraph &g) const; + + //MSpline3 *fAbsCoeffOzone; + //MSpline3 *fAbsCoeffAerosols; + + TGraph *fAbsCoeffOzone; + TGraph *fAbsCoeffAerosols; + +public: + MAtmosphere(const MCorsikaRunHeader &h) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0) + { + Init(h);//, "ozone.txt", "aerosols.txt"); + } + + MAtmosphere(const char *name1=0, const char *name2=0) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0) + { + if (name1) + InitOzone(name1); + if (name2) + InitAerosols(name2); + } + + ~MAtmosphere(); + + Float_t GetWavelengthMin() const; + Float_t GetWavelengthMax() const; + + Bool_t HasValidOzone() const; + Bool_t HasValidAerosol() const; + + Bool_t IsAllValid() const { return IsValid() && HasValidOzone() && HasValidAerosol(); } + + void PreCalcOzone(); + void PreCalcAerosol(); + Bool_t InitOzone(const TString name=""); + Bool_t InitAerosols(const TString name=""); + + void Init(const MCorsikaRunHeader &h, const char *name1=0, const char *name2=0) + { + MAtmRayleigh::Init(h); + + InitOzone(name1); + InitAerosols(name2); + } + + Double_t CalcOzoneAbsorption(Double_t h, Double_t wavelength, Double_t theta) const; + Double_t CalcAerosolAbsorption(Double_t h, Double_t wavelength, Double_t theta) const; + Double_t GetTransmission(const MPhotonData &ph) const; +}; + +#endif Index: trunk/Mars/msim/MSimAtmosphere.cc =================================================================== --- trunk/Mars/msim/MSimAtmosphere.cc (revision 19762) +++ trunk/Mars/msim/MSimAtmosphere.cc (revision 19763) @@ -53,8 +53,10 @@ #include "MPhotonData.h" +#include "MAtmosphere.h" + ClassImp(MSimAtmosphere); using namespace std; - +/* // ========================================================================== // @@ -253,41 +255,4 @@ } - /* - // The "orginal" code for the vertical thickness - Double_t GetVerticalThickness(Double_t obslev, Double_t height) const - { - // This is a C++-version of the original code from attenu.c - - // Limits (height in cm) of the four layers in which atmosphere is parametrized: - const double lahg[5] = - { - obslev, - TMath::Max(obslev, 4e5), - 1.0e6, - 4.0e6, - 1.0e7 - }; - - Double_t Rho_Tot = 0.0; - for (int i=0; i<4; i++) - { - const Double_t b = fAtmB[i]; - const Double_t c = fAtmC[i]; - - const Double_t h0 = TMath::Min(height, lahg[i+1]); - const Double_t h1 = lahg[i]; - - // Calculate rho for the i-th layer from the lower - // to the higher layer boundary. - // If height is within the layer only calculate up to height. - Rho_Tot += b*(exp(-h1/c) - exp(-h0/c)); - - if (lahg[i+1] > height) - break; - } - - return Rho_Tot; - } - */ Double_t CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const { @@ -334,5 +299,5 @@ // Calculate the traversed "vertical thickness" of air using the // US Standard atmosphere: - const Double_t Rho_tot = GetVerticalThickness(/*fObsLevel,*/ height); + const Double_t Rho_tot = GetVerticalThickness(height); // We now convert from "vertical thickness" to "slanted thickness" @@ -387,25 +352,4 @@ // checks are done for each call! return g.GetN()==0 ? 0 : g.Eval(wavelength)*1e-5; // from km^-1 to cm^-1 -/* - // Linear interpolation - // (FIXME: Is it faster to be replaced with a binary search?) - // ( This might be faster because we have more photons - // with smaller wavelengths) - //int index; - //for (index = 1; index =2 && i+1 km - rc += oz_conc[i] * 1e5/cos(theta); - - rc -= oz_conc[2]*0.2*1e5/cos(theta); - rc += oz_conc[i+1]*fmod(height/1e5,1)*1e5/cos(theta); - - return rc; - } - */ Double_t CalcOzoneAbsorption(Double_t h, Double_t wavelength, Double_t theta) const @@ -593,5 +516,5 @@ return 1; - //******* Ozone absorption ******* + // ******* Ozone absorption ******* if (h > 50.e5) h = 50.e5; @@ -617,5 +540,5 @@ return 1; - //******* Mie (aerosol) ******* + // ******* Mie (aerosol) ******* if (h > 30.e5) h = 30.e5; @@ -665,13 +588,13 @@ const Double_t h = TMath::Sqrt(H*H + d*d + 2*H*z) - R(); - //**** Rayleigh scattering: ***** + // **** Rayleigh scattering: ***** const Double_t T_Ray = CalcTransmission(h, wavelength, sin2); if (T_Ray<0) return 0; - //****** Ozone absorption: ****** + // ****** Ozone absorption: ****** const Double_t T_Oz = CalcOzoneAbsorption(h, wavelength, theta); - //******** Mie (aerosol) ******** + // ******** Mie (aerosol) ******** const Double_t T_Mie = CalcAerosolAbsorption(h, wavelength, theta); @@ -690,5 +613,5 @@ 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}; - +*/ // ========================================================================== Index: trunk/Mars/msim/Makefile =================================================================== --- trunk/Mars/msim/Makefile (revision 19762) +++ trunk/Mars/msim/Makefile (revision 19763) @@ -29,5 +29,6 @@ MSimAbsorption.cc \ MSimPointingPos.cc \ - MClonesArray.cc + MClonesArray.cc \ + MAtmosphere.cc ############################################################ Index: trunk/Mars/msim/SimLinkDef.h =================================================================== --- trunk/Mars/msim/SimLinkDef.h (revision 19762) +++ trunk/Mars/msim/SimLinkDef.h (revision 19763) @@ -20,3 +20,6 @@ #pragma link C++ class MClonesArray-; // - needed as custom streamer is already implemented +#pragma link C++ class MAtmosphere+; +#pragma link C++ class MAtmRayleigh+; + #endif