GeometryGenerator.cc

Go to the documentation of this file.

/*
  GeometryGeneratorKG

  this is useful for CONEX profiles (maybe also for GH profile)
  to set a shower axis (theta, phi)
*/

#include <utl/config.h>

#include "GeometryGenerator.h"

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>

#include <det/Detector.h>

#include <evt/Event.h>
#include <evt/Header.h>
#include <evt/ShowerSimData.h>

#include <sevt/SEvent.h>
#include <sevt/Header.h>

#include <fevt/FEvent.h>
#include <fevt/Header.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>

#include <utl/AugerCoordinateSystem.h>
#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/ErrorLogger.h>
#include <utl/RandomEngine.h>
#include <utl/Reader.h>
#include <utl/UTMPoint.h>
#include <utl/Vector.h>
#include <utl/Point.h>

#include <sdet/SDetector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>

#include <CLHEP/Random/RandFlat.h>

#include <math.h>
#include <cstddef>
#include <iostream>
#include <string>

using namespace utl;
using namespace fwk;
using namespace std;
using CLHEP::RandFlat;


namespace GeometryGeneratorKG {

  VModule::ResultFlag
  GeometryGenerator::Init()
  {
    CentralConfig& cc = *CentralConfig::GetInstance();

    Branch topBranch = cc.GetTopBranch("GeometryGenerator");

    // check for azimuth sampling
    // (CONEX profiles may be re-aligned in azimuth)
    Branch bAzimuth = topBranch.GetChild("SampleAzimuth");
    if (bAzimuth) {
      fSampleAzimuth = true;
      bAzimuth.GetChild("azimuthMin").GetData(fAzimuthMin);
      bAzimuth.GetChild("azimuthMax").GetData(fAzimuthMax);
    }

    // check for zenith sampling ------------------
    Branch bZenith = topBranch.GetChild ("SampleZenith");
    if (bZenith) {
      fSampleZenith = true;
      bZenith.GetChild("zenithMin").GetData(fZenithMin);
      bZenith.GetChild("zenithMax").GetData(fZenithMax);
      bZenith.GetChild("n_cos").GetData(fNcos);

      if (fZenithMin < 0 || fZenithMax < fZenithMin) {
        ostringstream err;
        err << "inconsistent zenith range: zenithMin = "
            << fZenithMin / degree << ", zenithMax = " << fZenithMax / degree;
        ERROR(err);
        return eFailure;
      }
      if (fNcos != 0 && (fZenithMin <= 90*degree) != (fZenithMax <= 90*degree)) { // XOR
        ostringstream err;
        err << "Zenith range can be only down or up going in n_cos > 0: n_cos = " << fNcos << ", "
               "zenithMin = " << fZenithMin / degree << ", zenithMax = " << fZenithMax / degree;
        ERROR(err);
        return eFailure;
      }
      if (fNcos != 0 && fNcos != 1 && fZenithMin > 90*degree) { // fNcos not 0 Nor 1 and upgoing
        ostringstream err;
        err << "n_cos must be 0 or 1 for up-going: n_cos = " << fNcos << ", "
               "zenithMin = " << fZenithMin / degree << ", zenithMax = " << fZenithMax / degree;
        ERROR(err);
        return eFailure;
      }
    }

    fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

    return eSuccess;
  }


  VModule::ResultFlag
  GeometryGenerator::Run(evt::Event& event)
  {
    evt::ShowerSimData& shower = event.GetSimShower();

    if (shower.HasDirection()) {
      ERROR("Nothing to do. Shower already has a direction. Check your ModuleSequence!");
      return eFailure;
    }

    ostringstream info;

    if (fSampleAzimuth) {
      const double azimuth = GetPhi();
      shower.SetGroundParticleCoordinateSystemAzimuth(azimuth);
      info << " set azimuth=" << azimuth/deg << " deg ";
    }

    if (fSampleZenith) {
      const double zenith = GetTheta();
      shower.SetGroundParticleCoordinateSystemZenith(zenith);
      info << " set zenith=" << zenith/deg << " deg ";
    }

    INFO(info);
    return eSuccess;
  }


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


  // GetTheta(): dice a zenith angle from isotropic flux on flat surface:
  // dN/dcos(theta) ~ cos(theta)

  double
  GeometryGenerator::GetTheta()
  {
    if (fZenithMin == fZenithMax)
      return fZenithMax;

    const auto n1 = fNcos + 1;
    const double t1 = pow(cos(fZenithMax), n1);
    const double t2 = pow(cos(fZenithMin), n1);

    const double c = RandFlat::shoot(&fRandomEngine->GetEngine(), t1, t2);

    double theta;
    if (fNcos == 0)
      theta = acos(c);
    else if (fNcos == 1) {
      if (fZenithMin > 90*degree)
        theta = acos(-sqrt(c));
      else
        theta = acos(sqrt(c));
    } else
      theta = acos(exp(log(c) / n1));

    return theta;
  }


  // GetPhi(): dice an azimuth angle

  double
  GeometryGenerator::GetPhi()
  {
    if (fAzimuthMax == fAzimuthMin)
      return fAzimuthMin;

    return RandFlat::shoot(&fRandomEngine->GetEngine(), fAzimuthMin, fAzimuthMax);
  }

}