ReadAerosolHL.cc

Go to the documentation of this file.

#include <sstream>

#include "ReadAerosolHL.h"
#include <utl/ErrorLogger.h>
#include <utl/TabulatedFunctionErrors.h>

#include <evt/Event.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>

#include <atm/AttenuationResult.h>
#include <atm/ScatteringResult.h>

#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/Point.h>

#include <atm/MonthlyAvgDBProfileModel.h>

using namespace std;
using namespace utl;
using namespace fwk;    
using namespace det;
using namespace fdet;
using namespace atm;

using namespace ReadAerosolHLNS;

ReadAerosolHL::ReadAerosolHL(){}

ReadAerosolHL::~ReadAerosolHL(){}

VModule::ResultFlag ReadAerosolHL::Init(){
  //
  //
  INFO("ReadAerosolHL::Init()");
  return eSuccess;
}

VModule::ResultFlag ReadAerosolHL::Run(evt::Event& event){
  //
  //
  INFO("ReadAerosolHL::Run()");

    Detector& det = Detector::GetInstance();

    // Set detector time to something
    //
    //    det.Update(TimeStamp(2005, 5, 11, 9, 0));
    det.Update(TimeStamp(799837213));

    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();

    // Define a track segment by specifying two points in the sky.  In this
    // example, we specify the points in the coordinate system of Los Leones
    //
    CoordinateSystemPtr 
      csLL(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetEyeCoordinateSystem() );
    Point A(100.*m,  1.*km,  1.*km,   csLL);
    Point B(2500.*m, 1. *km, 50.*m,   csLL);
    
    // compute Mie attenuation between points A and B.  For this calculation
    // we pass a set of wavelengths at which we want the computation to be carried out.
    //
    vector<double> wlength;
    wlength.push_back(300.*nanometer);
    wlength.push_back(350.*nanometer);
    wlength.push_back(400.*nanometer);

    AttenuationResult attRes;
    try {
      attRes = theAtm.EvaluateMieAttenuation(A, B, wlength);
    } catch (AugerException& ex) {
      cout << "An exception was thrown with the message " << ex.GetMessage() << endl;
      exit(EXIT_FAILURE);
    }

    cout << endl;
    cout << "Transmission coefficient at 300 nm = " << 
      attRes.GetTransmissionFactor().Y(300.*nanometer) << endl;

    cout << "Transmission coefficient at 350 nm = " << 
      attRes.GetTransmissionFactor().Y(350.*nanometer) << endl;
    
    cout << "Transmission coefficient at 375 nm = " << 
      attRes.GetTransmissionFactor().Y(375.*nanometer) << endl;

    cout << "Transmission coefficient at 400 nm = " << 
      attRes.GetTransmissionFactor().Y(400.*nanometer) << endl;

      
    // Compute Mie-scattered photon intensity at angle alpha,
    // and distance dist away from the track segment defined
    // by points A and B.

    ScatteringResult scatRes;
    try {
      scatRes = theAtm.EvaluateMieScattering(A, B, 30.*degree, 500.*m, wlength);
    }
    catch (...) {
      cout << "Could not do the computation. " << endl;
      exit(EXIT_FAILURE);
    }   
    cout << "scattered photon intensity at 350 nm = " << 
      scatRes.GetScatteringFactor().Y(350.*nanometer)*meter*meter 
         << " photons/m^2" << endl;      
    cout << endl;


    // Note that if you request data for a time when there is 
    // no data in the database, a NoDataForModelException will
    // be thrown.  You should check for this possibility. 
    //
    Detector::GetInstance().Update(UTCDateTime(2001, 1, 1).GetTimeStamp()); // no aerosol data for 2001 !
    
    try {
      theAtm.EvaluateMieAttenuation(A, B, wlength);
    }
    catch (NoDataForModelException& ex) {
      cout << "An exception was thrown when trying to retrieve" << endl;
      cout << "aerosol data for 2001.  This is probably because no data " << endl;
      cout << "exists in the database for this time." << endl;
      cout << "Note that if you change the MieModel type in AtmosphereInterfaceConfig " << endl;
      cout << "to 'ParametricXML' or 'Super' rather than 'MeasuredDB', then the " << endl;
      cout << "parametric model will be used, and you should NOT see this message" << endl;
    }   
    
    // If you request data for a portion of the sky for which
    // there is no valid data for the current time, you will 
    // get an out of bound exception
    //
    det.Update( TimeStamp(799837213));  // reset time to epoch when there is some data
    Point C(0.*m, 0.*m, 100.*m,   csLL);
    Point D(0.*m, 0.*m, 50.*km,   csLL);
    
    AttenuationResult wontWork;
    try {
      wontWork = theAtm.EvaluateMieAttenuation(C, D, wlength);
    } catch (AugerException& ex) {
      cout << "An expected exception of type " << ex.GetExceptionName() << " was thrown. " << endl;
      cout << "Message from the exception: " << ex.GetMessage() << endl;
    }


  return eSuccess;
}

VModule::ResultFlag ReadAerosolHL::Finish(){
  //
  //
  INFO("ReadAerosolHL::Finish()");
  return eSuccess;
}