Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 9350) +++ trunk/MagicSoft/Mars/Changelog (revision 9351) @@ -53,4 +53,18 @@ * mimage/MHHillas.cc, mimage/MMHillas.Ext.cc, mimage/MHHillasSrc.cc: - converted default binning from mm to deg + + * mraw/MRawRunHeader.cc: + - print less when values doen't have any meaning + + * msim/MSimAtmosphere.[h,cc]: + - added + + * msim/SimLinkDef.h, msim/Makefile: + - added MSimAtmosphere + + * mjobs/MJSimulation.cc: + - added simulation of atmosphere + - clean isolated pixels away before calculating image parameters + - don't clean calibration runs Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 9350) +++ trunk/MagicSoft/Mars/NEWS (revision 9351) @@ -75,5 +75,13 @@ memory and it was never freed, so the program took this memory until the end) - + + * Converted the code from MMCS' reflector program to MSimAtmosphere + (instead of the tables build in for MAGIC Winter and Summer + atmosphere we use the values written by Corsika fitted to the + tabulated atmosphere given as input) + + * For caluclating the displayed image parameters in ceres a basic + image cleaning removing isolated pixels is now done. + ;star Index: trunk/MagicSoft/Mars/mjobs/MJSimulation.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSimulation.cc (revision 9350) +++ trunk/MagicSoft/Mars/mjobs/MJSimulation.cc (revision 9351) @@ -80,4 +80,5 @@ #include "MSimMMCS.h" #include "MSimAbsorption.h" +#include "MSimAtmosphere.h" #include "MSimReflector.h" #include "MSimPointingPos.h" @@ -326,4 +327,5 @@ MSimMMCS simmmcs; + MSimAtmosphere simatm; MSimAbsorption absapd("PhotonDetectionEfficiency"); MSimAbsorption absmir("MirrorReflectivity"); @@ -511,6 +513,8 @@ // ------------------------------------------------------------------- - // FIXME: Remove isolated pixels - MImgCleanStd clean(7, 4.5); + // Remove isolated pixels + MImgCleanStd clean(0, 0); + clean.SetCleanLvl0(0); // The level above which isolated pixels are kept + clean.SetCleanRings(0); clean.SetMethod(MImgCleanStd::kAbsolute); @@ -551,5 +555,6 @@ // if (header.IsPointRun()) // tasks.AddToList(&stars); - //tasks.AddToList(&print); + if (1) + tasks.AddToList(&simatm); // Here because before fillh1 tasks.AddToList(&fillh1); tasks.AddToList(&fillG); @@ -558,4 +563,6 @@ tasks.AddToList(&absapd); tasks.AddToList(&absmir); + if (0) + tasks.AddToList(&simatm); // FASTER? } tasks.AddToList(&reflect); @@ -628,5 +635,6 @@ tasks.AddToList(&fillx0a); tasks.AddToList(&fillx0c); - //tasks.AddToList(&clean); + if (!header.IsCalibrationRun()) + tasks.AddToList(&clean); tasks.AddToList(&hcalc); tasks.AddToList(&fillx0d); Index: trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc =================================================================== --- trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc (revision 9350) +++ trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc (revision 9351) @@ -957,6 +957,9 @@ if (fFormatVersion>5) *fLog << "Observation: '" << fObservationMode << "'" << endl; - *fLog << "Source: '" << fSourceName << "' " << " "; - *fLog << fSourceEpochChar << dec << fSourceEpochDate << endl; + if (fSourceName[0]!=0 || fSourceEpochChar[0]!=0 || fSourceEpochDate!=0) + { + *fLog << "Source: '" << fSourceName << "' " << " "; + *fLog << fSourceEpochChar << dec << fSourceEpochDate << endl; + } if (fRunStart) *fLog << "Run Start: " << fRunStart << endl; @@ -974,8 +977,11 @@ else *fLog << (fNumSamplesLoGain+fNumSamplesHiGain) * fNumBytesPerSample << "B/pix" << endl; - *fLog << "Evt Counter: " << fNumEvents; - if (fFormatVersion>8) - *fLog << " (read=" << fNumEventsRead << ")"; - *fLog << endl; + if (fNumEvents>0 || fNumEventsRead>0) + { + *fLog << "Evt Counter: " << fNumEvents; + if (fFormatVersion>8) + *fLog << " (read=" << fNumEventsRead << ")"; + *fLog << endl; + } if (TString(t).Contains("header", TString::kIgnoreCase)) Index: trunk/MagicSoft/Mars/msim/MSimAtmosphere.cc =================================================================== --- trunk/MagicSoft/Mars/msim/MSimAtmosphere.cc (revision 9351) +++ trunk/MagicSoft/Mars/msim/MSimAtmosphere.cc (revision 9351) @@ -0,0 +1,869 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MSimAtmosphere.h" + +#include + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCorsikaRunHeader.h" +#include "MPhotonEvent.h" +#include "MPhotonData.h" + +ClassImp(MSimAtmosphere); + +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. +// +class MAtmRayleigh +{ +private: + static const Double_t fgMeanFreePath; // [g/cm^2] Mean free path for scattering Rayleigh XR + + 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[4]; // The index refers to the atmospheric layer (starting from sea level and going upwards) + Float_t fAtmB[4]; // The index refers to the atmospheric layer (starting from sea level and going upwards) + Float_t fAtmC[4]; // The index refers to the atmospheric layer (starting from sea level and going upwards) + + Double_t fRho[5]; // Precalculated integrals for rayleigh scatterning + + // -------------------------------------------------------------------------- + // + // 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 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 + + const Double_t h[5] = + { + fObsLevel, // 0km + TMath::Max(fObsLevel, 4e5), // 4km + 1.0e6, // 10km + 4.0e6, // 40km + 1.0e7 // 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; + } + } + +protected: + Double_t fObsLevel; // [cm] observation level a.s.l. + +public: + // Init an atmosphere from the data stored in MCorsikaRunHeader + MAtmRayleigh(const MCorsikaRunHeader &h) + { + Init(h); + } + + // Defualt constructor + MAtmRayleigh() : fObsLevel(-1) { } + + // 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; } + + // Init an atmosphere from the data stored in MCorsikaRunHeader + // This initialized fObsLevel, fR, fAtmB and fAtmC and + // PreCalcRho + void Init(const MCorsikaRunHeader &h) + { + // Observation level above earth radius + fObsLevel = h.GetObsLevel(); + + // Use earth radius as defined in Corsika + fR = h.EarthRadius(); + + //memcpy(fAtmA, (Float_t*)h.GetAtmosphericCoeffA(), sizeof(Float_t)*4); + memcpy(fAtmB, (Float_t*)h.GetAtmosphericCoeffB(), sizeof(Float_t)*4); + memcpy(fAtmC, (Float_t*)h.GetAtmosphericCoeffC(), sizeof(Float_t)*4); + + PreCalcRho(); + } + + // 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 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 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); + } +}; + +// ========================================================================== + +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 + { + // 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 GetWavelengthMax() const { return fAbsCoeffOzone && fAbsCoeffAerosols ? TMath::Min(fAbsCoeffOzone->GetX()[fAbsCoeffOzone->GetN()-1], fAbsCoeffAerosols->GetX()[fAbsCoeffAerosols->GetN()-1]) : -1; } + + Bool_t HasValidOzone() const { return fAbsCoeffOzone && fAbsCoeffOzone->GetN()>0; } + Bool_t HasValidAerosol() const { return fAbsCoeffAerosols && fAbsCoeffAerosols->GetN()>0; } + + Bool_t IsAllValid() const { return IsValid() && HasValidOzone() && HasValidAerosol(); } + + void 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 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 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 InitAerosols(const TString name="") + { + if (!name.IsNull()) + { + if (fAbsCoeffAerosols) + delete fAbsCoeffAerosols; + + fAbsCoeffAerosols = new TGraph(name); + fAbsCoeffAerosols->Sort(); + } + + if (!HasValidAerosol()) + return kFALSE; + + if (IsValid()) + PreCalcOzone(); + + return kTRUE; + } + + void Init(const MCorsikaRunHeader &h, const char *name1=0, const char *name2=0) + { + MAtmRayleigh::Init(h); + + InitOzone(name1); + InitAerosols(name2); + } +/* + 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 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 CalcAerosolAbsorption(Double_t h, Double_t wavelength, Double_t theta) const + { + if (!fAbsCoeffAerosols) + return 1; + + //******* Mie (aerosol) ******* + if (h > 30.e5) + h = 30.e5; + + // Mie absorption coefficient at sea level through interpolation + //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 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; + } +}; + +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}; + +// ========================================================================== + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MSimAtmosphere::MSimAtmosphere(const char* name, const char *title) + : fEvt(0), fAtmosphere(0) +{ + fName = name ? name : "MSimAtmosphere"; + fTitle = title ? title : "Simulate the wavelength and height-dependant atmpsheric absorption"; + + fAtmosphere = new MAtmosphere; +} + +// -------------------------------------------------------------------------- +// +// Calls Clear() +// +MSimAtmosphere::~MSimAtmosphere() +{ + delete fAtmosphere; +} + +// -------------------------------------------------------------------------- +// +// Search for the needed parameter containers. Read spline from file +// calling ReadFile(); +// +Int_t MSimAtmosphere::PreProcess(MParList *pList) +{ + fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent"); + if (!fEvt) + { + *fLog << err << "MPhotonEvent not found... aborting." << endl; + return kFALSE; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +Bool_t MSimAtmosphere::ReInit(MParList *pList) +{ + MCorsikaRunHeader *h = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); + if (!h) + { + *fLog << err << "MCorsikaRunHeader not found... aborting." << endl; + return kFALSE; + } + + //if (fRunHeader->Has(MCorsikaRunHeader::kRefraction)) + // *fLog << inf << "Atmospheric refraction already applied in Corsika... skipping our own." << endl; + + // FIXME: Check wavelength range + + /* + if (h->GetWavelengthMin()GetXmin()) + *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds lower bound of spline." << endl; + + if (h->GetWavelengthMax()>fSpline->GetXmax()) + *fLog << warn << "WARNING - Upper bound of wavelength bandwidth exceeds upper bound of spline." << endl; + */ + + fAtmosphere->Init(*h); + + if (!fAtmosphere->IsAllValid()) + { + *fLog << err << "ERROR - Something with the atmoshere's initialization went wrong!" << endl; + return kFALSE; + } + + if (h->GetWavelengthMin()GetWavelengthMin()) + *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds valid range of atmosphere." << endl; + + if (h->GetWavelengthMax()>fAtmosphere->GetWavelengthMax()) + *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds valid range of atmosphere." << endl; + + if (!h->Has(MCorsikaRunHeader::kAtmext)) + *fLog << warn << "WARNING - ATMEXT option not used for Corsika data." << endl; + + if (!h->Has(MCorsikaRunHeader::kRefraction)) + *fLog << warn << "WARNING - Refraction calculation disabled for Corsika data." << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +Int_t MSimAtmosphere::Process() +{ + // Get the number of photons in the list + const Int_t num = fEvt->GetNumPhotons(); + + // FIMXE: Add checks for + // * upgoing particles + // * Can we take the full length until the camera into account? + + // Counter for number of total and final events + Int_t cnt = 0; + for (Int_t i=0; iGetTransmission(ph); + + // Get a random value between 0 and 1 to determine whether the photon will survive + // gRandom->Rndm() = [0;1[ + if (gRandom->Rndm()>=eff) + continue; + + // Copy the surviving events bakc in the list + (*fEvt)[cnt++] = ph; + } + + // Now we shrink the array to the number of new entries. + fEvt->Shrink(cnt); + + return kTRUE; +} + +/* + Int_t MSimWavelength::Process() + { + // Get the number of photons in the list + const Int_t num = fEvt->GetNumPhotons(); + + // FIMXE: Add checks for + // * upgoing particles + // * wavelength range + // * check if corsika atmosphere is switched on + // * Can we take the full length until the camera into account? + + // Counter for number of total and final events + Int_t cnt = 0; + for (Int_t i=0; iGetWavelengthMin(); // WAVLGL + const Double_t max = fRunHeader->GetWavelengthMax(); // WAVLGU + const Double_t f = (max-min)/max; + + // WAVELENGTH = 1. / (1/min - RD(1)/(min*max/(max-min))) + + + ph.SetWavelength(TMath::Nint(min / (1. - gRandom->Rndm()*f))); + } + + return kTRUE; + } + */ + +// -------------------------------------------------------------------------- +// +// FileName: reflectivity.txt +// UseTheta: No +// +Int_t MSimAtmosphere::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; +/* + if (IsEnvDefined(env, prefix, "FileName", print)) + { + rc = kTRUE; + SetFileName(GetEnvValue(env, prefix, "FileName", fFileName)); + } + + if (IsEnvDefined(env, prefix, "UseTheta", print)) + { + rc = kTRUE; + SetUseTheta(GetEnvValue(env, prefix, "UseTheta", fUseTheta)); + } +*/ + return rc; +} Index: trunk/MagicSoft/Mars/msim/MSimAtmosphere.h =================================================================== --- trunk/MagicSoft/Mars/msim/MSimAtmosphere.h (revision 9351) +++ trunk/MagicSoft/Mars/msim/MSimAtmosphere.h (revision 9351) @@ -0,0 +1,38 @@ +#ifndef MARS_MSimAtmosphere +#define MARS_MSimAtmosphere + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MParList; +class MAtmosphere; +class MPhotonEvent; + +class MSimAtmosphere : public MTask +{ +private: + MPhotonEvent *fEvt; //! Event stroing the photons + + MAtmosphere *fAtmosphere; //! Instance of class describing atmosphere + + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + + // MTask + Int_t PreProcess(MParList *pList); + Bool_t ReInit(MParList *pList); + Int_t Process(); + +public: + MSimAtmosphere(const char *name=NULL, const char *title=NULL); + ~MSimAtmosphere(); + + // MSimAtmosphere + + // TObject + + ClassDef(MSimAtmosphere, 0) // Simulate the wavelength and height-dependant atmpsheric absorption +}; + +#endif Index: trunk/MagicSoft/Mars/msim/Makefile =================================================================== --- trunk/MagicSoft/Mars/msim/Makefile (revision 9350) +++ trunk/MagicSoft/Mars/msim/Makefile (revision 9351) @@ -26,4 +26,5 @@ MHPhotonEvent.cc \ MSimMMCS.cc \ + MSimAtmosphere.cc \ MSimAbsorption.cc \ MSimPointingPos.cc Index: trunk/MagicSoft/Mars/msim/SimLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/msim/SimLinkDef.h (revision 9350) +++ trunk/MagicSoft/Mars/msim/SimLinkDef.h (revision 9351) @@ -14,4 +14,5 @@ #pragma link C++ class MSimPointingPos+; #pragma link C++ class MSimAbsorption+; +#pragma link C++ class MSimAtmosphere+; #pragma link C++ class MSimMMCS+;