testInclinedProfileModel.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
#include <atm/InclinedAtmosphericProfile.h>

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/TabulatedFunction.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>

#include <cppunit/extensions/HelperMacros.h>
#include <tst/Verify.h>

#include <iostream>

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


class TestInclinedProfileModel : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(TestInclinedProfileModel);
  CPPUNIT_TEST_EXCEPTION(TestAccessHeightException,
                         InclinedAtmosphericProfile::InclinedAtmosphereModelException);
  CPPUNIT_TEST(TestAccessHeight);
  CPPUNIT_TEST(TestAccessDepth);
  CPPUNIT_TEST(TestVerticalGeometry);
  //CPPUNIT_TEST(TestAntiVerticalGeometry);
  //CPPUNIT_TEST(TestHorizontalGeometry);
  CPPUNIT_TEST_SUITE_END();

private:
  void InitTestProfile() { }

public:
  void
  setUp()
    override
  {
    CentralConfig::GetInstance(BOOTSTRAPFILE);
    Detector::GetInstance().Update(UTCDateTime(2005, 1, 1, 0, 0, 0).GetTimeStamp());
    ErrorLogger::GetInstance().SetVerbosity(Verbosity::eVerbose);
  }

  void tearDown() override { }

  void
  TestAccessHeightException()
  {
    Detector::GetInstance().GetAtmosphere().EvaluateHeightVsSlantDepth();
  }

  void
  TestAccessDepth()
  {
    const auto& cs = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    const auto& atm = Detector::GetInstance().GetAtmosphere();
    const Point origin(0,0,0, cs);
    const Vector dir(0,0,-1, cs);  // oposite of axis
    const double xStep = 20*g/cm2;
    atm.InitSlantProfileModel(origin, -dir, xStep);
    atm.EvaluateHeightVsSlantDepth();
  }

  void
  TestAccessHeight()  // DV: ??? this looks exactly the same as TestAccessDepth() ???
  {
    const auto& cs = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    const auto& atm = Detector::GetInstance().GetAtmosphere();
    const Point origin(0,0,0, cs);
    const Vector dir(0,0,-1, cs);  // oposite of axis
    const double xStep = 20*g/cm2;
    atm.InitSlantProfileModel(origin, -dir, xStep);
    atm.EvaluateHeightVsSlantDepth();
  }

  void
  TestInitWithProfile()
  {
    const auto& cs = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    InitTestProfile();
    // vetical geometry
    const Point core(0,0,0, cs);
    const Vector dir(0,0,-1, cs);
    const double xStep = 20*g/cm2;
    const auto& atm = Detector::GetInstance().GetAtmosphere();
    atm.InitSlantProfileModel(core, -dir, xStep);
  }

  void
  TestVerticalGeometry()
  {
    InitTestProfile();
    const auto& det = det::Detector::GetInstance();
    const auto& atm = det.GetAtmosphere();
    const auto& depthProfile = atm.EvaluateDepthVsHeight();
    //const auto& heightProfile = atm.EvaluateHeightVsDepth();

    // vetical geometry
    const auto& cs = det.GetReferenceCoordinateSystem();
    const Point core(0,0,0, cs);
    const Vector dir(0,0,-1, cs);
    const double xStep = 20*g/cm2;

    atm.InitSlantProfileModel(core, -dir, xStep);

    const auto& heightVsSlant = atm.EvaluateDistanceVsSlantDepth();

    const unsigned int nBins = 10; //fProfDepth.size()/4;
    const double minSlantDepth = heightVsSlant.MinX();
    const double maxSlantDepth = heightVsSlant.MaxX();
    const double dSlant = (maxSlantDepth - minSlantDepth) / (nBins - 1);
    cout << " -- " << minSlantDepth << ' '
         << maxSlantDepth << ' ' << nBins << '\n';
    // loop along shower and check mapping
    for (unsigned int i = 0; i < nBins; ++i) {
      const double atmSlantDepth = minSlantDepth + i * dSlant;
      const double height = heightVsSlant.Y(atmSlantDepth);
      const double atmDepth = depthProfile.Y(height);
      cout << i << ' '
           << atmSlantDepth/g*cm2 << ' '
           << atmDepth/g*cm2 << ' '
           << abs(atmSlantDepth - atmDepth)/g*cm2
           << endl;
      CPPUNIT_ASSERT(Verify<CloseTo>(atmDepth, atmSlantDepth, 1e-1*g/cm2));
    }
  }

/*
  NOT CLEAR HOW TO HANDLE THAT
  WHERE SHOULD THE SHOWER HAVE ITS FIRST INTERACTION ???

    void testAntiVerticalGeometry() {

        cout << " testAntiVerticalGeometry " << endl;

        initTestProfile();

        CoordinateSystemPtr referenceCS =
            det::Detector::GetInstance().GetReferenceCoordinateSystem();

        const atm::ProfileResult depthProfile
            = Detector::GetInstance().GetAtmosphere().
            EvaluateDepthVsHeight();

        const atm::ProfileResult heightProfile
            = Detector::GetInstance().GetAtmosphere().
            EvaluateHeightVsDepth();


        // vetical geometry
        Point core (0,0,0,referenceCS);
        Vector dir (0,0,1,referenceCS);

        const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();
        theAtm.InitSlantProfileModel (core, dir, fProfDepth);


        cout << " ok " << endl;

        const ProfileResult& slantVsHeight =
            theAtm.EvaluateSlantDepthVsHeight();

        unsigned int nBins = fProfDepth.size();
        double minHeight = slantVsHeight.MinX();
        double maxHeight = slantVsHeight.MaxX();
        cout << " -- " << minHeight
             << " " << maxHeight
             << " " << nBins
             << endl;


        // loop along shower and check mapping
        for (unsigned int i=0; i<nBins; i++) {

            double atmHeight = minHeight + (maxHeight-minHeight)*i/(nBins-1);
            //cout << atmHeight/m << " " << endl;
            //cout << " 1 " << endl;
            double atmDepth = depthProfile.Y (atmHeight);
            //cout << " 2 " << endl;
            double depth = slantVsHeight.Y (atmHeight);
            //cout << " 3 " << endl;

            cout << depth/g*cm*cm << " " << atmDepth/g*cm*cm
                 << " " << (depth/g*cm*cm == atmDepth/g*cm*cm)
                 << endl;
            CPPUNIT_ASSERT (Test<CloseTo>(depth/g*cm*cm,
                                          atmDepth/g*cm*cm,
                                          1.e-3));


        }

    }
*/

  void
  TestHorizontalGeometry()
  {
    InitTestProfile();
    auto& det = det::Detector::GetInstance();
    [[maybe_unused]] const auto& cs = det.GetReferenceCoordinateSystem();
    auto& atm = det.GetAtmosphere();

    const auto& depthProfile = atm.EvaluateDepthVsHeight();
    cout << "depth profile: X " << depthProfile.MinX()/meter << ' ' << depthProfile.MaxX()/meter << "\n"
            "               Y " << depthProfile.Y(depthProfile.MinX())/(g/cm2) << ' ' << depthProfile.Y(depthProfile.MaxX())/(g/cm2)
         << endl;

    const auto& heightProfile = atm.EvaluateHeightVsDepth();
    cout << "height profile: X " << heightProfile.MinX()/(g/cm2) << ' ' << heightProfile.MaxX()/(g/cm2) << "\n"
            "                Y " << heightProfile.Y(heightProfile.MinX())/meter << ' ' << heightProfile.Y(heightProfile.MaxX())/meter
         << endl;

    //Commented for the time being
    /*// horizontal geometry
    const Point core(0,0,0, referenceCS);
    const Vector dir(0,1,0, referenceCS);

    const Atmosphere& atm = Detector::GetInstance().GetAtmosphere();
    theAtm.InitSlantProfileModel(core, dir, 20*g/cm2);

    cout << " ok " << endl;

    const ProfileResult& slantVsHeight = atm.EvaluateSlantDepthVsHeight();

    const unsigned int nBins = 10; //fProfDepth.size();
    const double minHeight = slantVsHeight.MinX();
    const double maxHeight = slantVsHeight.MaxX();
    cout << " -- " << minHeight << ' '
         << maxHeight << ' ' << nBins << endl;

    // loop along shower and check mapping
    for (unsigned int i = 0; i < nBins; ++i) {

      const double atmHeight = minHeight + (maxHeight - minHeight) * i / (nBins - 1);
      const double atmDepth = depthProfile.Y(atmHeight);
      const double depth = slantVsHeight.Y(atmHeight);

      cout << depth/g*cm2 << ' ' << atmDepth/g*cm2 << ' '
           << (depth/g*cm2 == atmDepth/g*cm2) << endl;

      CPPUNIT_ASSERT(Verify<CloseTo>(depth/g*cm2, atmDepth/g*cm2, 1e-1));

    }*/
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(TestInclinedProfileModel);