MuonProductionDepthFinder.cc

Go to the documentation of this file.

#include "MuonProductionDepthFinder.h"
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/AtmosphereParameters.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/PMTCalibData.h>
#include <sevt/PMTRecData.h>
#include <sevt/StationRecData.h>
#include <sevt/Header.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/STimeVariance.h>

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

#include <utl/MathConstants.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/AxialVector.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/AugerUnits.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/PhysicalConstants.h>
#include <utl/UTMPoint.h>
#include <utl/TimeStamp.h>
#include <utl/Trace.h>
#include <utl/ErrorLogger.h>
#include <utl/TraceAlgorithm.h>
#include <utl/AugerException.h>
#include <utl/UTCDateTime.h>

#include <iostream>
#include <string>
#include <cmath>
#include <iterator>

using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace det;
using namespace sdet;
using namespace utl;
using namespace atm;
using namespace std;
using namespace MuonProductionDepthFinderGL;


double
ZToDepth(const double h,
         const vector<double>& hl,
         const vector<double>& al, const vector<double>& bl, const vector<double>& cl)
{
  for (unsigned int i = 0; i < 4; ++i)
    if (hl[i] <= h && h < hl[i+1])
      return al[i] + bl[i] * exp(-h / cl[i]);

  const double hTopAtmosphere = al[4] * cl[4] / bl[4];
  if (hl[4] <= h && h < hTopAtmosphere)
    return al[4] - bl[4] * (h / cl[4]);

  return 0;
}


double
DepthToZ(const double x,
         const vector<double>& hl,
         const vector<double>& al, const vector<double>& bl, const vector<double>& cl)
{
  vector<double> xl(5);
  for (unsigned int i = 0; i < 5; ++i)
    xl[i] = ZToDepth(hl[i], hl, al, bl, cl);

  for (unsigned int i = 0; i < 4; ++i)
    if (xl[i] >= x && x > xl[i+1])
      return -cl[i] * log((x - al[i]) / bl[i]);

  if (xl[4] >= x && x > 0)
    return cl[4] * (al[4] - x) / bl[4];

  return 0;
}


/*function extracted from AIRES code*/
double
ZToSlantDepth(const double zcground, const double theta, const double z,
              const vector<double>& hl,
              const vector<double>& al, const vector<double>& bl, const vector<double>& cl,
              const bool curvature)
{
  double zinjection = (al[4] * cl[4]) / bl[4];

  const double x0 = ZToDepth(zinjection, hl, al, bl, cl);

  const double cosTheta = cos(theta);

  const double rEarth = 6371315.*meter;  // 6371315<-CORSIKA, AIRES->6370949 (m)
  const double rEarth2 = 2 * rEarth;

  // The following parameters were set to optimize calculations whilst
  // ensuring that the error in Xslant is approximately less than
  // 0.1 g/cm2 (Xslant > 300 g/cm2).

  const double coplanemin = 0.996;
  const double deltax0 = 0.4;
  const double xquantum = 0.000000000002;

  // Checking input arguments.
  if (cosTheta <= 0)
    cout << "ERROR: you are using an wrong zenith angle (> 90 deg)" << endl;

  if (cosTheta < coplanemin && curvature) {

    double zv =
      sqrt(z*z + 2 * z * (rEarth + zcground) * cosTheta + (rEarth + zcground) * (rEarth + zcground)) - rEarth;
    const double xvertical = ZToDepth(zv, hl, al, bl, cl);
    const bool xsign = (x0 > xvertical);

    double xv;
    double xvmax;
    if (xsign) {
      xv = xvertical;
      xvmax = x0;
    } else {
      xv = x0;
      xvmax = xvertical;
    }

    const double cozen2 = cosTheta * cosTheta;
    const double xvrfact0 = deltax0 * cozen2;
    const double xvrfact = 1 + xvrfact0;
    const double xvqx = max(0.0001 * xvmax, xquantum);
    const double xvl0 = xvqx / xvrfact0;
    const double zgrsz2 = 4 * (1 - cozen2) * (zcground + rEarth) * (zcground + rEarth);

    double xpath = 0;
    while (xv < xvmax) {
      const double xvp = xv;

      if (xv > xvl0) {
        xv *= xvrfact;
      } else {
        xv += xvqx;
      }
      if (xv > xvmax)
        xv = xvmax;

      const double zvp = zv;
      zv = DepthToZ(xv, hl, al, bl, cl);

      const double sizprim2 = zgrsz2 /((zv + zvp + rEarth2) * (zv + zvp + rEarth2));
      xpath += (xv - xvp) / sqrt(1 - sizprim2);
    }

    return xsign ? -xpath : xpath;

  } else {

    // Quasi vertical axes are processed with a plane Earth model.
    const double zv = z * cosTheta + zcground;
    const double xvertical = ZToDepth(zv, hl, al, bl, cl);
    return (xvertical - x0) / cosTheta;

  }

  return 0;
}


double
FLogP0(double logz)
{
  if (logz < 3.2)
    logz = 3.2;
  if (logz > 6)
    logz = 6;

  return -0.6085 + logz * (1.955 + logz * (-0.3299 + logz * 0.0186));
}


double
EpsilonRz(const double r, const double z)
{
  const double lgz = log10(z);
  const double p0 = pow(10, FLogP0(lgz));
  return p0 / pow(r, 1.1758);
}


VModule::ResultFlag
MuonProductionDepthFinder::Init()
{
  INFO("MuonProductionDepthFinder::Init()");

  CentralConfig* const theCentralConfig = CentralConfig::GetInstance();
  const Branch& topBranch = theCentralConfig->GetTopBranch("MuonProductionDepthFinder");
  topBranch.GetChild("useConstRadiusCut").GetData(fUseConstRadiusCut);
  topBranch.GetChild("radiusCut").GetData(fRadiusCut);
  topBranch.GetChild("signalThreshold").GetData(fSignalThreshold);
  topBranch.GetChild("timeShift").GetData(fTimeShift);
  topBranch.GetChild("signalMin").GetData(fSignalMin);
  topBranch.GetChild("numberOfBins").GetData(fNXmuBins);
  topBranch.GetChild("xMin").GetData(fXmuMin);
  topBranch.GetChild("xMax").GetData(fXmuMax);
  topBranch.GetChild("usePionPathCorrection").GetData(fUsePionPathCorrection);

  fTblFuncRadiusCut =
    TabulatedFunction(
      topBranch.GetChild("tblEnergy").Get<vector<double> >(),
      topBranch.GetChild("tblRadiusCut").Get<vector<double> >()
    );

  return eSuccess;
}


VModule::ResultFlag
MuonProductionDepthFinder::Run(Event& event)
{
  INFO("MuonProductionDepthFinder::Run()");

  if (!(event.HasSEvent() && event.HasRecShower() &&
        event.GetRecShower().HasSRecShower()))
    return eContinueLoop;

  InitCoordinateSys(event, 0, 0, 0, 0);
  SetMuonProductionDepthHist(event);

  return eSuccess;
}


VModule::ResultFlag
MuonProductionDepthFinder::Finish()
{
  INFO("MuonProductionDepthFinder::Finish()");

  return eSuccess;
}


void
MuonProductionDepthFinder::InitCoordinateSys(const Event& event,
                                             const double dX, const double dY,
                                             const double dTheta, const double dPhi)
{
  INFO("MuonProductionDepthFinder::InitCoordSys()");

  // nothing to do if there is no SEvent
  const evt::ShowerSRecData& showersrec  = event.GetRecShower().GetSRecShower();
  const utl::Vector& axis = showersrec.GetAxis();
  const Point originalcore = showersrec.GetCorePosition();
  fCoreTime = showersrec.GetCoreTime();

  CoordinateSystemPtr originallocalCS = LocalCoordinateSystem::Create(originalcore);

  // moving the core
  const Vector shiftcore(dX, dY, 0, originallocalCS);
  fCore = originalcore + shiftcore;
  fCoreTime -= axis * shiftcore / kSpeedOfLight;
  fLocalCS = LocalCoordinateSystem::Create(fCore);

  // moving the angle
  fTheta = axis.GetTheta(originallocalCS) + dTheta;
  fPhi = axis.GetPhi(originallocalCS) + dPhi;

  // rotating to showerCS
  const utl::CoordinateSystemPtr alignedCS =
    CoordinateSystem::RotationZ(fPhi, fLocalCS);
  fShowerCS = CoordinateSystem::RotationY(fTheta, alignedCS);

  if (!fUseConstRadiusCut) {
    if (showersrec.GetEnergy() < fTblFuncRadiusCut.GetX(0))
      fRadiusCut = fTblFuncRadiusCut.GetY(0);
    else if (showersrec.GetEnergy() > fTblFuncRadiusCut.GetX(fTblFuncRadiusCut.GetNPoints() - 1))
      fRadiusCut = fTblFuncRadiusCut.GetY(fTblFuncRadiusCut.GetNPoints() - 1);
    else
      fRadiusCut = fTblFuncRadiusCut.Y(showersrec.GetEnergy());
  }
}


void
MuonProductionDepthFinder::SetMuonProductionDepthHist(Event& event)
{
  INFO("MuonProductionDepthFinder::SetMuonProductionDepthHist()");

  const ReferenceEllipsoid& wgs84 =
    ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  sevt::SEvent& sevent = event.GetSEvent();
  const det::Detector& detector = det::Detector::GetInstance();
  const sdet::SDetector& sDetector = detector.GetSDetector();

  evt::ShowerSRecData& showersrec = event.GetRecShower().GetSRecShower();

  if (!showersrec.HasMPDHistogram())
    showersrec.MakeMPDHistogram(fNXmuBins, fXmuMin, fXmuMax);

  const bool hasSimEvent = event.HasSimShower();
  if (hasSimEvent) {
    const evt::ShowerSimData& showerssim = event.GetSimShower();
    const evt::AtmosphereParameters& atmospherePars = showerssim.GetAtmosphereParameters();
    fAATM = atmospherePars.GetA();
    fBATM = atmospherePars.GetB();
    fCATM = atmospherePars.GetC();
    fHLAY = atmospherePars.GetLayerAltitudes();
  } else {
    const Atmosphere& atmo = Detector::GetInstance().GetAtmosphere();
    try {
      atmo.InitSlantProfileModel(fCore, showersrec.GetAxis(), 10*g/cm2);
    } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
      ERROR(e.GetMessage());
      return;
    }
    fSlantDepthProfile = atmo.EvaluateSlantDepthVsDistance();
  }

  const double cosZenith = cos(fTheta);
  const double sinZenith = sin(fTheta);

  //Loop over stations
  for (SEvent::StationIterator sIt = sevent.StationsBegin();
       sIt != sevent.StationsEnd(); ++sIt) {

    sevt::Station& station = *sIt;
    const sdet::Station& dStation = sDetector.GetStation(*sIt);
    const double fadcBinSize = dStation.GetFADCBinSize(); // this is already in ns
    //for correcting by the track length
    const double tankHeight = dStation.GetHeight();
    const double tankRadius = dStation.GetRadius();
    const double vemCorrection =
      cosZenith + 2 * tankHeight * sinZenith/(kPi * tankRadius);

    if (!station.IsCandidate())  // do not work with bad stations
      continue;

    bool useItInTheMPD = true;
    if (station.IsLowGainSaturation())  // do not work with saturated stations
      useItInTheMPD = false;

    const Point& sPosition = dStation.GetPosition();
    const Vector sRelPos = sPosition - fCore;

    const double delta = sRelPos.GetZ(fShowerCS);
    const double r = sRelPos.GetRho(fShowerCS);
    const double r2 = r * r;

    const TimeStamp fArrivalTimeOfShowerPlane =
      fCoreTime - TimeInterval(delta/kSpeedOfLight);
    const double fStartTimeReferredToArrivalOfShowerPlane =
      (station.GetRecData().GetSignalStartTime() - fArrivalTimeOfShowerPlane).GetInterval();

    const unsigned int startBin = station.GetRecData().GetSignalStartSlot();
    const unsigned int endBin = station.GetRecData().GetSignalEndSlot();

    double zground = wgs84.PointToLatitudeLongitudeHeight(sPosition).get<2>();

    // only tanks far from the core are used in the reconstruction
    if (r < fRadiusCut)
      useItInTheMPD = false;

    // minimal signal required in a tank for the MPD recontruction
    if (station.GetRecData().GetTotalSignal() < fSignalMin)
      useItInTheMPD = false;

    const int numberOfPMTs = distance(sIt->PMTsBegin(), sIt->PMTsEnd());

    vector<pair<double, double> > depthSignalVector;

    for (sevt::Station::PMTIterator pIt = sIt->PMTsBegin();
         pIt != sIt->PMTsEnd(); ++pIt) {
      depthSignalVector.clear();

      if (!pIt->HasRecData())
        continue;

      PMTRecData& pmtRec = pIt->GetRecData();
      /*
        if component "eDirectLightSubtracted" exist, it is used for the MPD,
        if not, the standard total component is used
      */
      if (!pmtRec.HasVEMTrace(StationConstants::eTotal))
        continue;
      sevt::StationConstants::SignalComponent sigComp = StationConstants::eTotal;
      if (pmtRec.HasVEMTrace(StationConstants::eDirectLightSubtracted))
        sigComp = StationConstants::eDirectLightSubtracted;

      const TraceD& vemtraceTot = pmtRec.GetVEMTrace(sigComp);

      double peak = 0;
      for (unsigned int i = startBin; i < endBin; ++i)
        peak = max(peak, vemtraceTot[i]);
      /*
        Signal threshold = 15% of the maximum of the trace
       */
      fSignalThreshold = peak * 0.15;

      for (unsigned int i = startBin; i < endBin; ++i) {

        const double ct =
          (fStartTimeReferredToArrivalOfShowerPlane + (i+0.5 - startBin) * fadcBinSize - fTimeShift) *
          kSpeedOfLight;

        if (ct <= 0)
          continue;

        double zmdelta = 0.5 * (r2 / ct - ct);
        double z = zmdelta + delta;
        /*
          Iterative estimation of epsilon(r, z - Delta).
          This converges quickly and a couple of iterations are sufficient.
        */
        for (unsigned int j = 0; j < 2; ++j) {
          const double epsilon = EpsilonRz(r, zmdelta);
          const double dcte = ct - epsilon * 0.5 * r2 / (sqrt(r2 + zmdelta * zmdelta));
          zmdelta = 0.5 * (r2 / dcte - dcte);
          z = zmdelta + delta;
          if (zmdelta < 0)
            break;
        }
        if (zmdelta < 0)
          continue;

        /*Correction for the pion path: It assumes that only the parent/last pion
          has transverse momentum (wrt the shower axis).
         */
        if (fUsePionPathCorrection) {
          const double energyPion = 1500.*GeV * pow(r, -0.84);  // L.Cazon Ph.D., ISBN 84-9750-467-4
          const double massPion = 139.57*MeV;
          const double ctauPion = 7.8045*meter;
          const double zPion = z / sqrt(r2 + z*z) * ctauPion * energyPion / massPion;
          z -= zPion;
        }

        if (vemtraceTot[i] > fSignalThreshold && z > 0) {

          if (hasSimEvent) {

            //For simulated events
            const bool curvature = false;
            fSlantDepth = ZToSlantDepth(zground, fTheta, z, fHLAY, fAATM, fBATM, fCATM, curvature);

          } else {

            //For real events
            const double atmMinDistance = fSlantDepthProfile.MinX() + 1*mm;
            const double atmMaxDistance = fSlantDepthProfile.MaxX() - 1*mm;
            const double distance = -z;
            if (distance < atmMinDistance || distance > atmMaxDistance)
              continue;
            fSlantDepth = fSlantDepthProfile.Y(distance);

          }

          const double signal = vemCorrection * vemtraceTot[i] / numberOfPMTs;
          if (useItInTheMPD)
            showersrec.GetMPDHistogram().Fill(fSlantDepth, signal);
          depthSignalVector.push_back(make_pair(fSlantDepth, signal));

        }
      } // loop over the trace

      if (!depthSignalVector.empty()) {
        if (!pmtRec.HasMuonProductionDepth())
          pmtRec.MakeMuonProductionDepth();
        TabulatedFunctionErrors& mpdProfile = pmtRec.GetMuonProductionDepth();
        for (unsigned int i = 0, n = depthSignalVector.size(); i < n; ++i) {
          // estimate of the uncertainty of the depth (preliminary)
          /*const double h0 = 7270*m;
          const double rho0 = 1.24e-3 * g/cm3;
          const double traceBinSize = sdet::Station::GetFADCBinSize();
          const double deltat = traceBinSize / sqrt(12.);
          double c = kSpeedOfLight;
          double x = depthVector[i];
          double deltaX = (x*h0*2/(cosZenith*r2)*pow(log(x*cosZenith/(h0*rho0)), 2)*c*deltat);*/
          const double deltaX = 0;
          mpdProfile.PushBack(depthSignalVector[i].first, deltaX, depthSignalVector[i].second, 0);
        }
      }
      pmtRec.SetIsUsedInGlobalMPD(useItInTheMPD);
    } // loop over the pmts

  } // loop over stations
}


bool
MuonProductionDepthFinder::IsContained(const SEvent& sevent)
  const
{
  const SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  for (SEvent::ConstStationIterator sIt = sevent.StationsBegin(),
       end = sevent.StationsEnd(); sIt != end; ++sIt) {
    const sevt::Station& station = *sIt;
    if (!station.IsCandidate())
      continue;
    const sdet::Station& dStation = sDetector.GetStation(*sIt);
    const unsigned int neighbours = dStation.GetCrown(1).size();
    if (neighbours < 6)
      return false;
  }
  return true;
}