AtmosphericProfile.cc

Go to the documentation of this file.

#include "AtmosphericProfile.h"

// basic headers with top-level offline classes
#include <evt/Event.h>
#include <fwk/CentralConfig.h>
#include <det/Detector.h>
#include <utl/ErrorLogger.h>

// geometry headers
#include <utl/UTMPoint.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <utl/CoordinateSystemPtr.h>

// These are the ones needed for the profiles
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

#include <TFile.h>
#include <TGraph.h>

using namespace std;
using namespace evt;
using namespace atm;
using namespace fwk;
using namespace utl;

using namespace AtmosphericProfileNS;

AtmosphericProfile::AtmosphericProfile(){}

AtmosphericProfile::~AtmosphericProfile(){}


VModule::ResultFlag
AtmosphericProfile::Init()
{
  //CentralConfig* theConfig =  CentralConfig::GetInstance();
  return eSuccess;
}


VModule::ResultFlag
AtmosphericProfile::Run(evt::Event& /*event*/)
{
  TFile file("profiles.root", "RECREATE");

  CoordinateSystemPtr pampaCS =
    CoordinateSystemRegistry::Get("PampaAmarilla");

  const Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
  // These do not need to be initialized because they are not slanted:
  // const ProfileResult& densityVsHeight = theAtm.EvaluateDensityVsHeight();
  const ProfileResult& depthVsHeight = theAtm.EvaluateDepthVsHeight();
  // const ProfileResult& heightVsDepth = theAtm.EvaluateHeightVsDepth();


  const Point corePosition(0,0,0,pampaCS);
  const UTMPoint utm_core(corePosition, ReferenceEllipsoid::eWGS84);
  const double z_ground = depthVsHeight.Y(utm_core.GetHeight());

  // Let's make a zenith angle scan and make a few plots:
  double azimuth = 0;
  TGraph* depthVsZenith = new TGraph();
  TGraph* depthVsZenith2 = new TGraph();
  for (double zenith = 0; zenith < 80*deg; zenith += 5*deg) {
    const CoordinateSystemPtr showerCs =
      CoordinateSystem::RotationY(zenith,
                                  CoordinateSystem::RotationZ(azimuth,
                                                              pampaCS));
    const Vector showerDirection(0,0,-1, showerCs);
    // Important! One has to call this method before using any slant profile
    // Note the minus sign (showerDirection points in the negative-z direction).
    theAtm.InitSlantProfileModel(corePosition, -showerDirection, 20*g/cm/cm);
    // const ProfileResult& slantProfile_DistanceVsDepth = theAtm.EvaluateDistanceVsSlantDepth();
    const ProfileResult& slantProfile_DepthVsDistance = theAtm.EvaluateSlantDepthVsDistance();
    const double x_ground = slantProfile_DepthVsDistance.Y(0);

    depthVsZenith->SetPoint(depthVsZenith->GetN(), zenith/deg, x_ground/(g/cm2));
    depthVsZenith2->SetPoint(depthVsZenith->GetN(), zenith/deg, z_ground/(g/cm2)/cos(zenith));
  }
  depthVsZenith->Write("depthVsZenithCurveAtm");
  depthVsZenith2->Write("depthVsZenithFlatAtm");

  const double maxh = 30*km;
  const double minh = 0*m;
  const unsigned int N = 200;
  const double dh = (maxh - minh)/N;
  TGraph* verticalDepthProfile = new TGraph();
  for (unsigned int i = 0; i != N; ++i) {
    const double h = minh+dh*i;
    const double depth = depthVsHeight.Y(h)/(g/cm2);
    verticalDepthProfile->SetPoint(i, h/km, depth);
  }
  verticalDepthProfile->Write("verticalDepthProfile");

  // Now let's plot the slant depth as a function of the distance to the ground for a particular zenith angle
  CoordinateSystemPtr showerCs = CoordinateSystem::RotationY(85*deg, CoordinateSystem::RotationZ(azimuth, pampaCS));
  Vector showerDirection(0,0,-1, showerCs);
  theAtm.InitSlantProfileModel(corePosition, -showerDirection, 20*g/cm/cm);
  const ProfileResult& slantProfile_DepthVsDistance = theAtm.EvaluateSlantDepthVsDistance();

  TGraph* slantDepthPlot = new TGraph();
  const double dmax = maxh/cos(85.*utl::deg);
  for (int i = 0; i != int(dmax/dh); ++i) {
    // SD people: note the minus sign!
    slantDepthPlot->SetPoint(i, i*dh/km, slantProfile_DepthVsDistance.Y(-i*dh)/g*cm2);
  }
  slantDepthPlot->Write("slantDepthProfile");

  file.Close();

  return eSuccess;
}


VModule::ResultFlag
AtmosphericProfile::Finish(){
  return eSuccess;
}