SdSimpleSim.cc

Go to the documentation of this file.

/*
     StationTrigger ErrorCode
     0  ok
     2  no-data, but station ok

     53 strange
     56 strange
*/

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

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

#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <sevt/Station.h>
#include <sevt/StationRecData.h>
#include <sevt/StationSimData.h>
#include <sevt/StationTriggerData.h>
#include <sevt/EventTrigger.h>

#include <utl/AxialVector.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Reader.h>
#include <utl/RandomEngine.h>
#include <utl/UTMPoint.h>

#include <det/Detector.h>
#include <sdet/SDetector.h>

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

#include "SdSimpleSim.h"
#include "MuonTimeModel.h"

#include <CLHEP/Random/RandPoisson.h>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RandGauss.h>
#include <CLHEP/Random/RandFlat.h>

#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TProfile.h>

#include <boost/range/iterator_range.hpp>

#include <iostream>

using CLHEP::RandGauss;
using CLHEP::RandFlat;
using CLHEP::RandPoisson;

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


const double SdSimpleSim::fNerlingA1a =  6.42522;
const double SdSimpleSim::fNerlingA1b = -1.53183;
const double SdSimpleSim::fNerlingA2a =  168.168;
const double SdSimpleSim::fNerlingA2b = -42.1368;

//
static double kMeterPerVEM = 1.2*meter;      // [m/VEM]
static double kTrackPerGEV = 5.25*meter/GeV; //(see E. Zas, Malargue Nov 04)
//static double MeanEnergy = 25.*MeV;         // e+- / gammas
// 240.MeV v. muon dE (1.2m) 0.240GeV/1.2m -> 5m/GeV


SdSimpleSim::SdSimpleSim() :
  fRandomEngine(RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::ePhysics)),
  fSdId(0),
  // from analysis of run sd_2004_10_16_01h51.root
  // per definition the T1 rate is about 100-110Hz
  fNoiseRatePerStation(1.0065346445791242e+02 * hertz)
{ }


VModule::ResultFlag
SdSimpleSim::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("SdSimpleSim");
  if (!topB) {
    ERROR("Could not find branch SdSimpleSim");
    return eFailure;
  }

  fSdId = 0;
  topB.GetChild("verbosity").GetData(fVerbosity);

  topB.GetChild("useT2LifeTime").GetData(fUseT2LifeTime);
  if (fUseT2LifeTime) {
    INFO("Using T2 life manager for SD station status.");
  }

  topB.GetChild("noise").GetData(fNoise);
  topB.GetChild("MuonProductionHeightFromProfile").GetData(fMuonProductionHeightFromProfile);
  topB.GetChild("n_sample").GetData(fNSample);
  topB.GetChild("LateralCut").GetData(fLateralCut);
  topB.GetChild("forceClosestTankToTrigger").GetData(fForceClosestTankToTrigger);
  topB.GetChild("onlyClosestTank").GetData(fOnlyClosestTank);
  topB.GetChild("planeShowerFront").GetData(fPlaneFront);
  topB.GetChild("MuonRScale").GetData(fMuonRScale);
  topB.GetChild("em_age_fact").GetData(fEMAgeFactor);

  // --- optional stuff ----------------------------

  Branch TMSB =  topB.GetChild("TimeModelStat");
  if (TMSB) {
    INFO("Write time model statistics to root file.");

    string FileName;
    TMSB.GetData(FileName);

    TFile* const StatFile = TFile::Open(FileName.c_str(), "recreate");
    StatFile->cd();

    fLorenzo_vs_X = 0;
    fdNdX = 0;
    fCharged = 0;
    fMuon = 0;
    fElectron = 0;
    fGamma = 0;
    fdEdX = 0;
    fLorenzo = 0;
    fdNdz = 0;
    fCharged_vs_z = 0;
    fMuon_vs_z = 0;
    fElectron_vs_z = 0;
    fGamma_vs_z = 0;
    fdEdX_vs_z = 0;

    fStatTree = new TTree("stat", "time model statistics");
    fStatTree->Branch("Lorenzo_vs_X", "TProfile", &fLorenzo_vs_X, 32000, 0);
    fStatTree->Branch("dNdX", "TProfile", &fdNdX, 32000, 0);
    fStatTree->Branch("Charged", "TProfile", &fCharged, 32000, 0);
    fStatTree->Branch("Gamma", "TProfile", &fGamma, 32000, 0);
    fStatTree->Branch("Electron", &fElectron, 32000, 0);
    fStatTree->Branch("Muon", "TProfile", &fMuon, 32000, 0);
    fStatTree->Branch("dEdX", "TProfile", &fdEdX, 32000, 0);
    fStatTree->Branch("Lorenzo", "TProfile", &fLorenzo, 32000, 0);
    fStatTree->Branch("dNdz", "TProfile", &fdNdz, 32000, 0);
    fStatTree->Branch("Charged_vs_z", "TProfile", &fCharged_vs_z, 32000, 0);
    fStatTree->Branch("Gamma_vs_z", "TProfile", &fGamma_vs_z, 32000, 0);
    fStatTree->Branch("Electron_vs_z", &fElectron_vs_z, 32000, 0);
    fStatTree->Branch("Muon_vs_z", "TProfile", &fMuon_vs_z, 32000, 0);
    fStatTree->Branch("dEdX_vs_z", "TProfile",  &fdEdX_vs_z, 32000, 0);

  } else {

    fStatTree = 0;

  }

  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "            verbosity: " << fVerbosity << "\n"
          "     use T2 lift time: " << fUseT2LifeTime << "\n"
          "       generate noise: " << fNoise << "\n"
          "   plane shower front: " << fPlaneFront << "\n"
          "   muon prod. profile: " << fMuonProductionHeightFromProfile << "\n"
          "    max no. samplings: " << fNSample << "\n"
          "          lateral cut: " << fLateralCut/km << "km\n"
          "   force closest tank: " << fForceClosestTankToTrigger << "\n"
          "   only closeset tank: " << fOnlyClosestTank << "\n"
          "    muon radial scale: " << fMuonRScale/m << "m\n"
          "      e.m. age factor: " << fEMAgeFactor << "\n";
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
SdSimpleSim::Run(evt::Event& event)
{
  // check event structure
  if (!event.HasSimShower())
    return eContinueLoop;

  const evt::ShowerSimData& SimShower = event.GetSimShower();

  if (!SimShower.HasDirection())
    return eContinueLoop;

  if (!SimShower.HasPosition())
    return eContinueLoop;

  if (!SimShower.HasLongitudinalProfile())
    return eContinueLoop;

  if (!SimShower.HasGHParameters())
    return eContinueLoop;

  const bool hasMuonProfile =
    SimShower.HasLongitudinalProfile(ShowerSimData::eMuon);

  const bool hasMuonProductionProfile =
    SimShower.HasLongitudinalProfile(ShowerSimData::eMuonProduction);

  const bool hasGammaProfile =
    SimShower.HasLongitudinalProfile(ShowerSimData::ePhoton);

  const bool hasElectronProfile =
    SimShower.HasLongitudinalProfile(ShowerSimData::eElectron);

  const bool hasdEdXProfile =
    SimShower.HasLongitudinalProfile(ShowerSimData::eEnergyDeposit);

  if (!hasMuonProfile || !hasGammaProfile || !hasElectronProfile) {
    ERROR("Input file not suitable for SdSimpleSim");
    return eContinueLoop;
  }

  if (!hasMuonProductionProfile && fMuonProductionHeightFromProfile) {
    ERROR("Cannot use MuonProductionHeight from simulation");
    return eContinueLoop;
  }

  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  // extrapolate CONEX profiles
  // assume exponential decay after the CONEX profile ends !
  // fit the decay after Xmax
  FitParam fitMuon = FitDecay(SimShower.GetLongitudinalProfile(ShowerSimData::eMuon));
  FitParam fitMuonProd;
  if (hasMuonProductionProfile)
    fitMuonProd = FitDecay(SimShower.GetLongitudinalProfile(ShowerSimData::eMuonProduction));
  FitParam fitGamma = FitDecay(SimShower.GetLongitudinalProfile(ShowerSimData::ePhoton));
  FitParam fitElectron = FitDecay(SimShower.GetLongitudinalProfile(ShowerSimData::eElectron));
  FitParam fitdEdX = FitDecay(SimShower.GetLongitudinalProfile(ShowerSimData::eEnergyDeposit), true);
  fProfileExtrapolation[ShowerSimData::eMuon] = fitMuon;
  fProfileExtrapolation[ShowerSimData::eMuonProduction] = fitMuonProd;
  fProfileExtrapolation[ShowerSimData::ePhoton] = fitGamma;
  fProfileExtrapolation[ShowerSimData::eElectron] = fitElectron;
  fProfileExtrapolation[ShowerSimData::eEnergyDeposit] = fitdEdX;
  // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  // make new time
  const utl::TimeStamp& coreTime = SimShower.GetTimeStamp();
  utl::TimeStamp TriggerTime = coreTime;
  bool First = true;

  // get/create the sd-event
  if (!event.HasSEvent())
    event.MakeSEvent();
  sevt::SEvent& sEvent = event.GetSEvent();
  sevt::Header& sEventHeader = sEvent.GetHeader();
  sEventHeader.SetTime(coreTime);
  sEventHeader.SetId(++fSdId);

  // get the detector definitions
  det::Detector& Det = det::Detector::GetInstance();
  const sdet::SDetector& SDet = Det.GetSDetector();
  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

  //const ProfileResult& heightVsSlantDepth = Det.GetAtmosphere().EvaluateHeightVsSlantDepth();
  //const ProfileResult& distanceVsHeight = Det.GetAtmosphere().EvaluateDistanceVsHeight();

  // get the simulated shower info
  double Xmax = SimShower.GetGHParameters().GetXMax();
  float simEnergy = SimShower.GetEnergy();
  int simPrimary = SimShower.GetPrimaryParticle();

  if (fVerbosity > 3) {
    ostringstream dbg;
    dbg << " Sim-Parameters:  Xmax= "<<Xmax/g*cm2;
    INFO(dbg);
  }

  const utl::Vector& simAxis = SimShower.GetDirection();
  const utl::Point& simCore = SimShower.GetPosition();
  const utl::TabulatedFunction& SimProfile = SimShower.GetLongitudinalProfile();
  const utl::CoordinateSystemPtr& ShowerCS = SimShower.GetShowerCoordinateSystem();

  const CoordinateSystemPtr localCS = SimShower.GetLocalCoordinateSystem();
  const double Zenith = (-SimShower.GetDirection()).GetTheta(localCS); 

  const double cosTheta = cos(Zenith);
  const double absCosTheta = std::fabs(cosTheta);
  //double sinTheta = sin(Zenith);

  // ------- if inclined atmosphere was not initialized before, do it here ---------
  try {
    Det.GetAtmosphere().EvaluateSlantDepthVsDistance();
  } catch (const InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
    Det.GetAtmosphere().InitSlantProfileModel(simCore, -simAxis, 5.*g/cm2);
  }
  // -------------------------------------------------------------------------------

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

  const atm::ProfileResult temperatureProfile = Det.GetAtmosphere().EvaluateTemperatureVsHeight();

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

  //  const ProfileResult& slantDepthVsHeight =
  const ProfileResult& slantDepthVsDistance = Det.GetAtmosphere().EvaluateSlantDepthVsDistance();

  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  // interpolate Atmosphere
  // assume exponential increase of slant depth with decreasing height
  fSlantDepthExtrapolation = FitAtmosphere(slantDepthVsDistance);
  // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  /*
    double showerXstart = (SimShower.GetFluorescencePhotons(0))[0].X();

    // find the part, that was emitting (fluorescence) light
    double InitialHeight = heightProfile.Y(showerXstart*cosTheta);

    const utl::Point showerStart = simCore +
    simAxis * InitialHeight/simAxis.GetZ(siteCS);
  */

  const utl::TabulatedFunction& SimMuonProfile =
    SimShower.GetLongitudinalProfile(ShowerSimData::eMuon);

  const utl::TabulatedFunction& SimGammaProfile =
    SimShower.GetLongitudinalProfile(ShowerSimData::ePhoton);

  const utl::TabulatedFunction SimElectronProfile =
    SimShower.GetLongitudinalProfile(ShowerSimData::eElectron);

  utl::TabulatedFunction SimMuonProductionProfile;
  if (hasMuonProductionProfile && fMuonProductionHeightFromProfile)
    SimMuonProductionProfile = SimShower.GetLongitudinalProfile(ShowerSimData::eMuonProduction);
  else {
    if (!fPlaneFront)
      ERROR("No muon production profile  !!!!!!!!!");
  }

  // just for debug plottings -------------------
  utl::TabulatedFunction SimdEdXProfile;
  if (hasdEdXProfile)
    SimdEdXProfile =
      SimShower.GetLongitudinalProfile(ShowerSimData::eEnergyDeposit);
  // --------------------------------------------

  // Maximum profile slant-depth
  double ProfileMaxDepth = SimProfile[SimProfile. GetNPoints() - 1].X();

  // ---------INFO-----------
  if (fVerbosity > 3) {
    ostringstream dbg;
    dbg << " max prof depth " << ProfileMaxDepth/g*cm*cm
        << " " << SimMuonProfile[0].X()/g*cm2
        << " " << SimMuonProfile[SimMuonProfile.GetNPoints()-1].X()/g*cm2;
    INFO(dbg);
  }
  // --------------------

  // calculate
  const double HeightObs = wgs84.PointToLatitudeLongitudeHeight(simCore).get<2>();

  // atmospheric parameters
  const double atmDistanceMin = slantDepthVsDistance.MinX();
  const double atmDistanceMax = slantDepthVsDistance.MaxX();
  const double AtmDepthMin = depthProfile.Y(depthProfile.MaxX());
  const double AtmDepthMax = depthProfile.Y(depthProfile.MinX());
  const double XobsVert = depthProfile.MinX() < HeightObs ?
    depthProfile.Y(HeightObs) :
    AtmDepthMax;
  const double Temperature =
    temperatureProfile.Y(fmax(depthProfile.MinX(), HeightObs));
  static double g_earth = 9.81 * m/(s*s);
  const double Pressure = XobsVert * g_earth;

  int lorenzoPrimary = -1;
  if (!fPlaneFront) {
    switch (simPrimary) {
    case utl::Particle::ePhoton:
      lorenzoPrimary = 0;
      break;
    case utl::Particle::eProton:
      lorenzoPrimary = 1;
      break;
    case utl::Particle::eIron:
      lorenzoPrimary = 2;
      break;
    default:
      ERROR("Unknown lorenzoPrimary!");
      return eFailure;
      break;
    }
  }

  // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  if (fStatTree) {

    INFO(" STAT TREE ");

    fLorenzo_vs_X = 0;
    fdNdX = 0;
    fCharged = 0;
    fMuon = 0;
    fElectron = 0;
    fGamma = 0;
    fdEdX = 0;
    fLorenzo = 0;
    fdNdz = 0;
    fCharged_vs_z = 0;
    fMuon_vs_z = 0;
    fElectron_vs_z = 0;
    fGamma_vs_z = 0;
    fdEdX_vs_z = 0;

    int nBin;
    double xMin, xMax;
    double zMin, zMax;
    double Xvert, HeightX, z;

    // charged
    nBin = SimProfile.GetNPoints();
    xMin = SimProfile[0].X();
    xMax = SimProfile[nBin-1].X();
    fCharged = new TProfile("charged", "charged profile",
                            nBin, xMin/g*cm*cm, xMax/g*cm*cm);
    Xvert = xMin * absCosTheta;
    if (Xvert < AtmDepthMin)
      Xvert = AtmDepthMin;
    if (Xvert > AtmDepthMax)
      Xvert = AtmDepthMax;
    HeightX = heightProfile.Y(Xvert);
    zMax = (HeightX - HeightObs) / absCosTheta;
    Xvert = xMax*absCosTheta;
    if (Xvert < AtmDepthMin)
      Xvert = AtmDepthMin;
    if (Xvert > AtmDepthMax)
      Xvert = AtmDepthMax;
    HeightX = heightProfile.Y(Xvert);
    zMin = (HeightX - HeightObs) / absCosTheta;
    fCharged_vs_z = new TProfile("charged_vs_z", "charged profile vs z",
                                 nBin, zMin/meter, zMax/meter);

    for (TabulatedFunction::ConstIterator iBin = SimProfile.Begin();
         iBin != SimProfile.End(); ++iBin) {

      fCharged->Fill(iBin->X()/g*cm2, iBin->Y());

      Xvert = iBin->X() * absCosTheta;
      if (Xvert < AtmDepthMin || Xvert > AtmDepthMax)
        continue;
      HeightX = heightProfile.Y(Xvert);
      z = (HeightX - HeightObs) / absCosTheta;

      fCharged_vs_z->Fill(z/meter, iBin->Y());

    }

    // gammas
    if (hasGammaProfile) {
      nBin = SimGammaProfile.GetNPoints();
      xMin = SimGammaProfile[0].X();
      xMax = SimGammaProfile[nBin-1].X();
      fGamma = new TProfile("gamma", "gamma profile",
                            nBin, xMin/g*cm*cm, xMax/g*cm*cm);
      Xvert = xMin * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMax = (HeightX - HeightObs) / absCosTheta;
      Xvert = xMax * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMin = (HeightX - HeightObs) / absCosTheta;
      fGamma_vs_z = new TProfile("gamma_vs_z", "gamma profile vs z",
                                 nBin, zMin/meter, zMax/meter);

      for (utl::TabulatedFunction::ConstIterator x = SimGammaProfile.Begin();
           x != SimGammaProfile.End(); ++x) {

        fGamma->Fill(x->X()/g*cm2, x->Y());

        Xvert = x->X() * absCosTheta;
        if (Xvert < AtmDepthMin || Xvert > AtmDepthMax)
          continue;
        HeightX = heightProfile.Y(Xvert);
        z = (HeightX - HeightObs) / absCosTheta;

        fGamma_vs_z->Fill(z/meter, x->Y());

      }
    }


    // electrons
    if (hasElectronProfile) {
      nBin = SimElectronProfile.GetNPoints();
      xMin = SimElectronProfile[0].X();
      xMax = SimElectronProfile[nBin-1].X();
      fElectron = new TProfile("electron", "electron profile",
                               nBin, xMin/g*cm*cm, xMax/g*cm*cm);
      Xvert = xMin * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMax = (HeightX - HeightObs) / absCosTheta;
      Xvert = xMax * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMin = (HeightX - HeightObs) / absCosTheta;
      fElectron_vs_z = new TProfile("electron_vs_z",
                                    "electron profile vs z",
                                    nBin, zMin/meter, zMax/meter);

      for (utl::TabulatedFunction::ConstIterator x = SimElectronProfile.Begin();
           x != SimElectronProfile.End(); ++x) {

        fElectron->Fill(x->X()/g*cm2, x->Y());

        Xvert = x->X() * absCosTheta;
        if (Xvert < AtmDepthMin || Xvert > AtmDepthMax)
          continue;
        HeightX = heightProfile.Y(Xvert);
        z = (HeightX - HeightObs) / absCosTheta;

        fElectron_vs_z->Fill(z/meter, x->Y());

      }
    }

    // muons
    if (hasMuonProfile) {
      nBin = SimMuonProfile.GetNPoints();
      xMin = SimMuonProfile[0].X();
      xMax = SimMuonProfile[nBin-1].X();
      fMuon = new TProfile("muon", "muon profile",
                           nBin, xMin/g*cm*cm, xMax/g*cm*cm);
      Xvert = xMin * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMax = (HeightX - HeightObs) / absCosTheta;
      Xvert = xMax * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMin = (HeightX - HeightObs) / absCosTheta;
      fMuon_vs_z = new TProfile("muon_vs_z", "muon profile vs z",
                                nBin, zMin/meter, zMax/meter);

      for (utl::TabulatedFunction::ConstIterator x = SimMuonProfile.Begin();
           x != SimMuonProfile.End(); ++x) {

        fMuon->Fill(x->X()/g*cm2, x->Y());

        Xvert = x->X() * absCosTheta;
        if (Xvert < AtmDepthMin || Xvert > AtmDepthMax)
          continue;
        HeightX = heightProfile.Y(Xvert);
        z = (HeightX - HeightObs) / absCosTheta;

        fMuon_vs_z->Fill(z/meter, x->Y());

      }
    }

    // dEdXs
    if (hasdEdXProfile) {
      nBin = SimdEdXProfile.GetNPoints();
      xMin = SimdEdXProfile[0].X();
      xMax = SimdEdXProfile[nBin-1].X();
      fdEdX = new TProfile("dEdX", "dEdX profile", nBin, xMin/g*cm2, xMax/g*cm2);
      Xvert = xMin * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMax = (HeightX - HeightObs) / absCosTheta;
      Xvert = xMax * absCosTheta;
      if (Xvert < AtmDepthMin)
        Xvert = AtmDepthMin;
      if (Xvert > AtmDepthMax)
        Xvert = AtmDepthMax;
      HeightX = heightProfile.Y(Xvert);
      zMin = (HeightX - HeightObs) / absCosTheta;
      fdEdX_vs_z = new TProfile("dEdX_vs_z", "dEdX profile vs z",
                                nBin, zMin/meter, zMax/meter);

      for (utl::TabulatedFunction::Iterator x = SimdEdXProfile.Begin();
           x != SimdEdXProfile.End(); ++x) {

        fdEdX->Fill(x->X()/g*cm2, x->Y());

        Xvert = x->X() * absCosTheta;
        if (Xvert < AtmDepthMin || Xvert > AtmDepthMax)
          continue;
        HeightX = heightProfile.Y(Xvert);
        z = (HeightX - HeightObs) / absCosTheta;

        fdEdX_vs_z->Fill(z/meter, x->Y());

      }
    }

#if 0
    // Lorenzos TimeModel   -------------------------
    if (hasMuonProductionProfile) {
      /*TimeModel tmpModel(Zenith/deg, SimMuonProductionProfile,
                           XobsVert, ShowerSimData::eMuonProduction);
      fdNdz = tmpModel.Get_dNdz_FromProfile();
      fdNdX = tmpModel.Get_dNdX_FromProfile(absCosTheta, HeightObs);*/
    } else if (hasMuonProfile) {
      /*TimeModel tmpModel(Zenith/deg, SimMuonProfile, XobsVert,
                           ShowerSimData::eMuon);
      fdNdz = tmpModel.Get_dNdz_FromProfile();
      fdNdX = tmpModel.Get_dNdX_FromProfile(absCosTheta, HeightObs);*/
    }

    /*TimeModel tmpModel(Zenith/deg, lorenzoPrimary);
    fLorenzo = tmpModel.Get_dNdz(zMin, zMax, nBin);
    fLorenzo_vs_X = tmpModel.Get_dNdX(zMin, zMax, nBin, absCosTheta, HeightObs);*/
#endif

    if (fLorenzo_vs_X == 0)
      fLorenzo_vs_X = new TProfile();

    if (fdNdX == 0)
      fdNdX = new TProfile();

    if (fCharged == 0)
      fCharged = new TProfile();

    if (fMuon == 0)
      fMuon = new TProfile();

    if (fElectron == 0)
      fElectron = new TProfile();

    if (fGamma == 0)
      fGamma = new TProfile();

    if (fdEdX == 0)
      fdEdX = new TProfile();

    if (fLorenzo == 0)
      fLorenzo = new TProfile();

    if (fdNdz == 0)
      fdNdz = new TProfile();

    if (fCharged_vs_z == 0)
      fCharged_vs_z = new TProfile();

    if (fMuon_vs_z == 0)
      fMuon_vs_z = new TProfile();

    if (fElectron_vs_z == 0)
      fElectron_vs_z = new TProfile();

    if (fGamma_vs_z == 0)
      fGamma_vs_z = new TProfile();

    if (fdEdX_vs_z == 0)
      fdEdX_vs_z = new TProfile();

    fStatTree->Fill();

    delete fLorenzo_vs_X;
    delete fdNdX;
    delete fCharged;
    delete fMuon;
    delete fElectron;
    delete fGamma;
    delete fdEdX;
    delete fLorenzo;
    delete fdNdz;
    delete fCharged_vs_z;
    delete fMuon_vs_z;
    delete fElectron_vs_z;
    delete fGamma_vs_z;
    delete fdEdX_vs_z;

  }
  // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  // Lorenzos TimeModel
  MuonTimeModel* lorenzoTimeModel = nullptr;
  if (!fPlaneFront) {

    if (fMuonProductionHeightFromProfile) {

      if (hasMuonProductionProfile) {
        
        if (fVerbosity > 2)
          INFO("TimeModel from MuonProduction Profile");
        TabulatedFunction logZdist =
          GetLogZdist(SimMuonProductionProfile, simCore, simAxis);

        lorenzoTimeModel = new MuonTimeModel(fRandomEngine, Zenith, &logZdist, true, true);

      } else if (hasMuonProfile) {

        WARNING ("TimeModel from Muon-Profile: NOT PROPERLY IMPLEMENTED YET");

        TabulatedFunction logZdist =
          CalculateLogZdist(SimMuonProfile, absCosTheta, XobsVert);

        lorenzoTimeModel = new MuonTimeModel(fRandomEngine, Zenith, &logZdist, true, true);

      } else {

        ERROR("no muon profiles available, to init time model!");
        return eFailure;

      }

    } else {
      if (fVerbosity > 2)
        INFO("TimeModel from parametrisation");
      const double logE = log10(simEnergy/eV);
      lorenzoTimeModel = new MuonTimeModel(fRandomEngine, Zenith, logE, lorenzoPrimary);
    }

    if (fVerbosity > 2)
      lorenzoTimeModel->Info();

  } // not plane shower front

  // calcualte S1000 ////////////////////

  const double XobsS1000 = XobsVert / absCosTheta;

  // if the profile range is exceeded, try fitting the range after maximum
  // if fit fails (for muon-number-profile) use simple muon decay of last entry

  double Ne = 0;
  double Ng = 0;
  double Nmu = 0;

  if (XobsS1000 > ProfileMaxDepth) {
    WARNING("Shower does not reach ground! Extrapolate profile!");

    if (fProfileExtrapolation[ShowerSimData::eMuon].p0 == 0 &&
        fProfileExtrapolation[ShowerSimData::eMuon].p1 == 0) { //extrapolation failed
      const double NmuLastPoint = SimMuonProfile[SimMuonProfile.GetNPoints()-1].Y();
      const double depthLastPoint = SimMuonProfile[SimMuonProfile.GetNPoints()-1].X()/g*cm2;
      const double distanceLastPoint = (log(depthLastPoint) - fSlantDepthExtrapolation.p0)/fSlantDepthExtrapolation.p1/m;
      const double gammaFactor = 10;

      Nmu = NmuLastPoint * exp((distanceLastPoint/m - HeightObs/m) / (kSpeedOfLight/m*s * kMuonLifetime/s * gammaFactor));

      if (fVerbosity > 3) {
        ostringstream dbg;
        dbg <<" Nmu-profile (last point)= "<<NmuLastPoint
            <<" Depth (last point)= "<<depthLastPoint
            <<" distance (last point)= "<<distanceLastPoint
            <<" height needed= "<<HeightObs
            <<" gammaFactor= "<<gammaFactor
            <<" --> Nmu= "<<Nmu;
        INFO(dbg);
      }
    } else {
      Nmu = exp(fProfileExtrapolation[ShowerSimData::eMuon].p0 +
                fProfileExtrapolation[ShowerSimData::eMuon].p1 * XobsS1000/g*cm2);
    }

    Ng = exp(fProfileExtrapolation[ShowerSimData::ePhoton].p0 +
             fProfileExtrapolation [ShowerSimData::ePhoton].p1 * XobsS1000/g*cm2);
    Ne = exp(fProfileExtrapolation[ShowerSimData::eElectron].p0 +
             fProfileExtrapolation [ShowerSimData::eElectron].p1 * XobsS1000/g*cm2);
  } else {
    Nmu = SimMuonProfile.Y(XobsS1000);
    Ng = SimGammaProfile.Y(XobsS1000);
    Ne = SimElectronProfile.Y(XobsS1000);
  }

  const double Age = utl::ShowerAge(XobsS1000, Xmax);
  const double Rm = utl::MoliereRadius(Temperature, Pressure, absCosTheta);

  // ---------INFO-----------
  if (fVerbosity > 3) {
    ostringstream dbg;
    dbg << " ********* calculate S1000 ********** " << endl;
    dbg << " Age " << Age << " rm " << Rm/m
        << " Ne " << Ne << " Ng " << Ng
        << " Nmu " << Nmu << " XobsS1000 " << XobsS1000/g*cm*cm;
    INFO(dbg);
  }
  // --------------------

  const double S1000 = CalculateTankSignal(1000.*m, 1.8*m, 1.2*m, Zenith,
                                           Age, Rm,
                                           Ne, Ng, Nmu, false);

  // ---------INFO-----------
  if (fVerbosity > 3) {
    ostringstream dbg;
    dbg << " S1000 " << S1000 << " Energy/EeV " << simEnergy/EeV << endl;
    dbg << " ********* calculated S1000 ********** ";
    INFO(dbg);
  }
  // --------------------

  int closestStationID = -1;

  // ---------------------------------------------------------------------
  // loop stations to find closest one
  if (fForceClosestTankToTrigger || fOnlyClosestTank) {

    double closestStationDistance = 0;

    for (sdet::SDetector::StationIterator iStation = SDet.StationsBegin();
         iStation != SDet.StationsEnd(); ++iStation) {

      const int tankID = iStation->GetId();

      const Point& tankPos = iStation->GetPosition();

      const double PosXinShower = tankPos.GetX(ShowerCS);
      const double PosYinShower = tankPos.GetY(ShowerCS);
      //double PosZinShower = tankPos.GetZ(ShowerCS);

      const double coreDistance = sqrt(PosXinShower*PosXinShower +
                                       PosYinShower*PosYinShower);

      if (closestStationID < 0 || coreDistance < closestStationDistance) {

        // skip stations not T2-alive
        if (fUseT2LifeTime) {
          if (iStation->IsInAcquisition()) {
            closestStationID = tankID;
            closestStationDistance = coreDistance;
          }
        } else {
          closestStationID = tankID;
          closestStationDistance = coreDistance;
        }
      }

    } // end loop stations
  }

  int nNotInDAQ = 0;

  // ----------------------------------------------------------------------
  // loop stations (detector definitions)
  for (sdet::SDetector::StationIterator iStation = SDet.StationsBegin();
       iStation != SDet.StationsEnd(); ++iStation) {

    // skip stations, and don't push it into the SEvent, if is was not T2-alive
    if (fUseT2LifeTime && !iStation->IsInAcquisition()) {
      ++nNotInDAQ;
      continue;
    }

    const int tankID = iStation->GetId();

    // ------------------------------------------------------------------
    //   Create station - event
    //   (ALL stations are part of the event !)
    //
    if (sEvent.HasStation(tankID)) {
      ERROR("StationID occured twice in detector !!!!!! skipping !!!");
      continue;
    }
    sEvent.MakeStation(tankID);
    sevt::Station& tank = sEvent.GetStation(tankID);

    // don't forget the tank trigger data
    // trigger is eNone and errorCode=1
    tank.MakeTriggerData();
    sevt::StationTriggerData& stationTrigger = tank.GetTriggerData();

//#warning How to properly calculate FD-T3 time offsets for the stations? Need FD T3!

    sevt::StationTriggerData::Algorithm tankTriggerAlgo = sevt::StationTriggerData::eNone;
    int errorCode = 2;
    int timeOffset = 0;
    double timeWindow = GetCDASTriggerTimeWindow(coreTime);
    stationTrigger.SetOffsetMicroSecond(timeOffset / microsecond);
    stationTrigger.SetWindowMicroSecond(timeWindow / microsecond);
    stationTrigger.SetAlgorithm(tankTriggerAlgo);
    stationTrigger.SetErrorCode(errorCode);
    //stationTrigger.SetPLDTrigger(const int trigger);

    // DV: setting stations to silent prevents later setting it to candidate;
    // all stations should be kept candidate and at the end all stations
    // wothout an signal should be set tank.SetSilent(); // all station are created 'silent'
    // ----------------------------------------------------------------

    // check if only closest station should be simulated
    if (fOnlyClosestTank && tankID != closestStationID)
      continue;

    const Point& tankPos = iStation->GetPosition();

    const double PosXinShower = tankPos.GetX(ShowerCS);
    const double PosYinShower = tankPos.GetY(ShowerCS);
    const double PosZinShower = tankPos.GetZ(ShowerCS);

    const double coreDistance = sqrt(PosXinShower*PosXinShower + PosYinShower*PosYinShower);

    /*if (coreDistance > fLateralCut) {
      continue;
    }*/

    const Point tankAtAxis(0,0, PosZinShower, ShowerCS);

    const double tankOverSeaLevel =
      wgs84.PointToLatitudeLongitudeHeight(tankPos).get<2>();

    const double tankAtAxisOverSeaLevel =
      wgs84.PointToLatitudeLongitudeHeight(tankAtAxis).get<2>();

    const Vector tankAtAxisToCore = tankAtAxis - simCore;
    const double distanceTankAtAxisToCore = tankAtAxisToCore.GetMag();
    const double distanceFromCoreAtAxis = distanceTankAtAxisToCore *
      (tankAtAxisToCore*simAxis < 0 ? 1 : -1);

    //double tankOverSeaLevel = UTMPoint(tankPos, wgs84).GetHeight();
    //double ThetaOftankInShower = (showerStart - tankPos).GetTheta(siteCS);

    // ---------INFO-----------
    if (fVerbosity > 3) {
      ostringstream dbg;
      dbg  << " tank " << tankID
           << " height@axis " << tankAtAxisOverSeaLevel
           << " AtmMin " << atmDistanceMin
           << " AtmMax " << atmDistanceMax;
      INFO(dbg);
    }
    // --------------------

    bool flagAtmExceeded = false;
    double xTank = 0;
    if (distanceFromCoreAtAxis < atmDistanceMin ||
        distanceFromCoreAtAxis > atmDistanceMax) {

      // ---------INFO-----------
      if (fVerbosity > 3) {
        ostringstream dbg;
        dbg  << " ATM EXCEEDED: "
             << " TankR: " << coreDistance/m
             << " TankAtAxisOverSeaLevel: " << tankAtAxisOverSeaLevel/m
             << " DistanceAtAxisFromCore: " << distanceFromCoreAtAxis/m
             << " min " << atmDistanceMin << " " << atmDistanceMax;
        INFO(dbg);
      }
      // --------------------

      xTank = exp(fSlantDepthExtrapolation.p0 + fSlantDepthExtrapolation.p1 * distanceFromCoreAtAxis) * g/cm2;

      flagAtmExceeded = true;

    } else {

      //tankAtAxisOverSeaLevel);
      xTank = slantDepthVsDistance.Y(distanceFromCoreAtAxis);

    }

    const double tankRadius = iStation->GetRadius();
    const double tankHeight = iStation->GetHeight();

    if (coreDistance > fLateralCut) {
      //tank.SetSilent();
      continue;
    }

    double Xobs = xTank;

    bool ShowerProfileRangeExceeded = false;
    if (Xobs > ProfileMaxDepth) {

      ShowerProfileRangeExceeded = true;
      //Xobs = ProfileMaxDepth;
    }

    // ---------INFO-----------
    if (fVerbosity > 1) {
      ostringstream dbg;
      dbg  << " ShowerProfileRangeExceeded " << ShowerProfileRangeExceeded
           << " maxX " << ProfileMaxDepth/g*cm*cm
           << " Xobs " << Xobs/g*cm*cm
           << " Xobs " << Xobs
           << " " << SimMuonProfile[0].X()
           << " " << SimMuonProfile[SimMuonProfile.GetNPoints()-1].X()
           << " " << SimMuonProfile[0].X()/g*cm2
           << " " << SimMuonProfile[SimMuonProfile.GetNPoints()-1].X()/g*cm2;
      INFO(dbg);
    }
    // --------------------

    // ---------INFO-----------
    if (fVerbosity > 1 && ShowerProfileRangeExceeded) {
      ostringstream dbg;
      dbg  << " ele "
           << " p0 " << fProfileExtrapolation [ShowerSimData::eElectron].p0
           << " p1 " << fProfileExtrapolation [ShowerSimData::eElectron].p1
           << " gam "
           << " p0 " << fProfileExtrapolation [ShowerSimData::ePhoton].p0
           << " p1 " << fProfileExtrapolation [ShowerSimData::ePhoton].p1
           << " muo "
           << " p0 " << fProfileExtrapolation [ShowerSimData::eMuon].p0
           << " p1 " << fProfileExtrapolation [ShowerSimData::eMuon].p1;
      INFO(dbg);
    }
    // --------------------

    // if the profile range is exceeded, try fitting the range after maximum
    // if fit fails (for muon-number-profile) use simple muon decay of last entry
    double Ne = 0;
    double Ng = 0;
    double Nmu= 0;
    if (ShowerProfileRangeExceeded) {
      if (fProfileExtrapolation[ShowerSimData::eMuon].p0 == 0 &&
          fProfileExtrapolation[ShowerSimData::eMuon].p1 == 0) { //extrapolation failed
        const double NmuLastPoint = SimMuonProfile[SimMuonProfile.GetNPoints()-1].Y();
        const double depthLastPoint = SimMuonProfile[SimMuonProfile.GetNPoints()-1].X()/g*cm2;
        const double distanceLastPoint = (log(depthLastPoint) - fSlantDepthExtrapolation.p0)/fSlantDepthExtrapolation.p1/m;
        const double gammaFactor = 10;

        Nmu = NmuLastPoint *
          exp((distanceLastPoint/m - distanceFromCoreAtAxis/m)/(kSpeedOfLight/m*s * kMuonLifetime/s *gammaFactor));

        if (fVerbosity > 3) {
          ostringstream dbg;
          dbg << " Nmu-profile (last point)= " << NmuLastPoint
              << " Depth (last point)= " << depthLastPoint
              << " distance (last point)= " << distanceLastPoint
              << " distance needed at axis= " << distanceFromCoreAtAxis
              << " gammaFactor= " << gammaFactor
              << " --> Nmu= " << Nmu;
          INFO(dbg);
        }
      } else {
        Nmu = exp(fProfileExtrapolation[ShowerSimData::eMuon].p0 +
                  fProfileExtrapolation[ShowerSimData::eMuon].p1 * Xobs/g*cm2);
      }

      Ng = exp(fProfileExtrapolation[ShowerSimData::ePhoton].p0 +
               fProfileExtrapolation[ShowerSimData::ePhoton].p1 * Xobs/g*cm2);
      Ne = exp(fProfileExtrapolation[ShowerSimData::eElectron].p0 +
               fProfileExtrapolation[ShowerSimData::eElectron].p1 * Xobs/g*cm2);
    } else {
      Nmu = SimMuonProfile.Y(Xobs);
      Ng = SimGammaProfile.Y(Xobs);
      Ne = SimElectronProfile.Y(Xobs);
    }

    const double Age = utl::ShowerAge(Xobs, Xmax);
    const double Rm = utl::MoliereRadius(Temperature, Pressure, absCosTheta);
    double tankSignal =
      CalculateTankSignal(coreDistance, tankRadius, tankHeight, Zenith, Age, Rm, Ne, Ng, Nmu);

    // ---------INFO-----------
    if (fVerbosity > 3) {
      ostringstream dbg;
      dbg  << " Age " << Age << " rm " << Rm/m
           << " Ne " << Ne << " Ng " << Ng
           << " Nmu " << Nmu << " Xtank " << Xobs/g*cm2
           << " Xmax " << Xmax/g*cm2
           << " S " << tankSignal;
      INFO(dbg);
    }
    // --------------------

    //bool tankIsForced = false; // unused

    if (tankSignal < 0.01) {
      if (fVerbosity > 3) {
        ostringstream dbg;
        dbg << " tank " << tankID
            << " has signal < 0.01 (closestStation= " << closestStationID << ")  ";
        INFO(dbg);
      }

      if (fForceClosestTankToTrigger && tankID == closestStationID) {
        if (fVerbosity > 1) {
          ostringstream infoString;
          infoString << "Forcing tank trigger for tank signal (previous) S=" << tankSignal << " set to S=5";
          //         tankIsForced = true; // unused
          INFO(infoString);
        }
        // use a fixed tank signal for non-triggered tanks with
        // fForceClosestTankToTrigger
        tankSignal = 5;
      } else {
        if (fVerbosity > 3) {
          INFO(" setting to silent");
        }
        //tank.SetSilent(); <-- should not be done here
        continue;
      }
    }

    // Shower plane - simulated arrival time

    double tankPlaneSimTime = - PosZinShower / kSpeedOfLight;

    // Set simulation truth for all tanks: time, signal (components)
    if (tank.HasSimData()) {
      ostringstream msg;
      msg << "Station " << tank.GetId() << " already has SimData!";
      WARNING(msg);
    } else {
      tank.MakeSimData();
    }
    sevt::StationSimData& tankSim = tank.GetSimData();
    tankSim.SetSimulatorSignature("SdSimpleSimKG");

    TimeInterval tankSimTime = tankPlaneSimTime;
    /*
    Particle muons(Particle::eMuon, Particle::eShower, tankPos, simAxis, tankSimTime, fNMuons, fSMuons);
    Particle electrons(Particle::eElectron, Particle::eShower, tankPos, simAxis, tankSimTime, fNElectrons, fSElectrons);
    Particle gammas(Particle::ePhoton, Particle::eShower, tankPos, simAxis, tankSimTime, fNGammas, fSGammas);
    tankSim.AddParticle(muons);
    tankSim.AddParticle(electrons);
    tankSim.AddParticle(gammas);
    */
    tankSim.SetPlaneFrontTime(coreTime + tankSimTime);


    const double p_trig = T1TriggerProbability(tankSignal, S1000, Zenith/deg);

    // ----------INFO-----------
    if (fVerbosity > 3)  {
      ostringstream dbg;
      dbg << " p_trig " << p_trig;
      INFO(dbg);
    }
    // ---------------------

    if (fForceClosestTankToTrigger && tankID == closestStationID) {
      ; // do nothing
    } else {
      if (RandFlat::shoot(&fRandomEngine.GetEngine(), 0.0, 1.0) > p_trig) {

        if (fVerbosity > 0) {
          ostringstream info;
          info << "tank id="<< tankID
               << " r=" << coreDistance/m << "m"
                  " S=" << tankSignal << "VEM"
                  " p_trig=" << p_trig
               << " (NO TRIGGER)";
          INFO(info);
        }

        //tank.SetSilent(); <-- should not be done here
        continue;
      }
    }

    //if (tankSignal>10)
    //tankTriggerAlgo = sevt::StationTriggerData::eThreshold;

    //if (ShowerProfileRangeExceeded || flagAtmExceeded) {
    if (flagAtmExceeded && fVerbosity > 2) {

      // if we really enter here, there is a station with a possible signal, that
      // can not get calculated from the given shower profile.

      ostringstream warn;
      warn << "Shower profile range exceeded "
              "(" << (ShowerProfileRangeExceeded ? "true" : "false")
           << ")! "
           << (flagAtmExceeded ? "[ATM out of bound]" : "")
           << " Station at X=" << xTank/g*cm2 << " g/cm2 "
              "(h=" << tankOverSeaLevel/m << "m,"
              "r=" << coreDistance/m << "m) "
              "but Profile ends at " << ProfileMaxDepth/g*cm2
           << "g/cm2 !!!\n"
           << "(Signal: " << tankSignal << ", "
              "p_trig: " << p_trig << ", "
              "nMuon: " << fNMuons << ", "
              "nElectrons: " << fNElectrons << ", "
              "nGammas: " << fNGammas << ")";
      WARNING(warn);

      //return eContinueLoop;
    }

    // Shower plane - approximation

    double tankTime = tankSimTime.GetInterval();

    /*
      from Lorenzo Cazon
    */
    const double R = tankPos.GetRho(ShowerCS);
    double Phi = tankPos.GetPhi(ShowerCS);
    if (!fPlaneFront) {

      while (Phi < 0)
        Phi += 360.*deg;
      lorenzoTimeModel->SetCoordinates(R/m, Phi/rad);

      //if (!tankIsForced) { // what is this good for ????
      if (fNMuons) {
//#warning Need confirmation from Lorenzo !!

        // limit number of muons
        int nDice = fNMuons;
        if (nDice > 1000)
          nDice = 1000;
        tankTime = lorenzoTimeModel->GetFirstTime(fNMuons);

      } else {

//#warning Need electron arrival time model here !!
        WARNING(" TIME-MODEL: Muon-model timing for signal without muons !!! Use nMuon=1");
        //continue;
        tankTime = lorenzoTimeModel->GetFirstTime(1.0);

      }
      //}
    }

    // Detector time-resolution
    // dependent on FADC trace (shower width, particle distribution,
    //                          FADC traces)
    // from GAP 2005-035
    static const double constK = 5.368e-5*ns/(m*m);

    static const double SigmaT0 = 26.84*ns;

    double sigmaTime = 0;
    if (!fPlaneFront)    // shower to shower fluctuations already from Lorenzo
      sigmaTime = SigmaT0;  // --> only detector resolution needs to be added
    else
      sigmaTime = sqrt(SigmaT0*SigmaT0 + pow(constK*R*R*absCosTheta, 2));

    double deltaT = RandGauss::shoot(&fRandomEngine.GetEngine(), 0, sigmaTime);

    tankTime += deltaT;
    // --- end of calculation of arrival times -----------

    utl::TimeStamp stationTime = coreTime + TimeInterval(tankTime);

    if (First) {
      First = false;
      TriggerTime = stationTime;
    } else {
      if (TriggerTime > stationTime)
        TriggerTime = stationTime;
    }

    tankTriggerAlgo = sevt::StationTriggerData::eTimeOverThreshold;
    errorCode = 0;
    stationTrigger.SetAlgorithm(tankTriggerAlgo);
    stationTrigger.SetErrorCode(errorCode);

    /*if (tankSignal < 1) { closestStationID
      tank.SetSilent(); // this is not done here, but in SdEventSelector
      continue;
    }*/

    // What about saturated stations?
    tank.SetHighGainSaturation(false);
    tank.SetLowGainSaturation(false);
    tank.SetCandidate();

    // If not silent, fill the rec-data
    tank.MakeRecData();
    sevt::StationRecData& tankRec = tank.GetRecData();

    tankRec.SetTotalSignal(tankSignal);
    tankRec.SetSignalStartTime(stationTime);

    if (fVerbosity > 0) {
      ostringstream info;
      info << "tank id="<< tankID
           << " r=" << coreDistance/m << "m"
           << " S=" << tankSignal << "VEM"
           << " p_trig=" << p_trig
           << " t=" << tankTime << "ns";
      if (tankID == closestStationID)
        info << " [closest]";
      if (fVerbosity > 3)
        info << " z_sh=" << PosZinShower/m << "m"
             << " z_sh/c= " << PosZinShower/kSpeedOfLight << "ns";
      INFO(info.str());

      // cout << endl;
    }

  } // loop stations

  delete lorenzoTimeModel;

  if (nNotInDAQ) {
    ostringstream info;
    info << "There have been " << nNotInDAQ << " stations not in DAQ accoring to T2 life file!";
    INFO(info);
  }

  if (!sEvent.HasTrigger())
    sEvent.MakeTrigger();

  sevt::EventTrigger& STrigger = sEvent.GetTrigger();
  STrigger.SetTime(TriggerTime);
  STrigger.SetAlgorithm("FD");
  //STrigger.SetSender("FD");
  //STrigger.SetId("FD");
  //STrigger.SetPreviousId("FD");
  //STrigger.SetNumberOfStations("FD");
  //STrigger.SetSDPAngle();

  // THE NOISE loop
  // loop all stations again
  if (fNoise) {

    // ---------INFO--------------
    if (fVerbosity > 3)
      INFO(" -------->noise ");
    // -----------------------

    int nNoiseTank = 0;
    int nSignalNoiseTank = 0;

    // -------------------------------------------------------------------------
    // noise loop
    for (sdet::SDetector::StationIterator iStation = SDet.StationsBegin();
         iStation != SDet.StationsEnd(); ++iStation) {

      // skip stations, which are not T2-alive
      if (fUseT2LifeTime) {
        if (!iStation->IsInAcquisition()) {
          continue;
        }
      }

      int tankID = iStation->GetId();

      // central trigger (station) time window
      double timeWindow = GetCDASTriggerTimeWindow(coreTime);
      double timeOffset = 0;

      // station is not part of event yet (-> cannot happen ... )
      // if station had no shower signal yet
      if (!sEvent.HasStation (tankID)) {

        if (Noise(sEvent, TriggerTime, tankID, timeOffset, timeWindow))
          ++nNoiseTank;

        continue;
      }

      sevt::Station& tank = sEvent.GetStation(tankID);

      // check for TRIGGER and adjust TimeWindow, TimeOffset
      if (true == false &&                                        // turn off noise for signal stations
          tank.HasTriggerData() &&
          (tank.GetTriggerData().GetAlgorithm() == sevt::StationTriggerData::eT1Threshold ||
           tank.GetTriggerData().GetAlgorithm() == sevt::StationTriggerData::eT2Threshold ||
           tank.GetTriggerData().GetAlgorithm() == sevt::StationTriggerData::eTimeOverThreshold)) {

        const double fadcTraceLength = iStation->GetFADCTraceLength() * iStation->GetFADCBinSize();
        timeWindow = fadcTraceLength;
        timeOffset = 0;
        if (Noise(sEvent, TriggerTime, tankID, timeOffset, timeWindow)) {
          ++nNoiseTank;
          ++nSignalNoiseTank;
        }

        continue;
      }

      // tank belong to event, but had no trigger yet
      if (Noise(sEvent, TriggerTime, tankID, timeOffset, timeWindow))
        ++nNoiseTank;

    } // end noise loop

    if (fVerbosity > 1) {
      ostringstream inf;
      inf << " Number of noise stations: " << nNoiseTank;
      if (nSignalNoiseTank > 0) {
        inf << ", of which " << nSignalNoiseTank << " tank" << (nSignalNoiseTank == 1 ? "" : "s")
            << " also contain" << " shower signal.";
      }
      INFO(inf);
    }

  } // end noise loop ----------------------------------------------------------

  // ---------INFO--------------
  if (fVerbosity > 3)
    INFO( "noise ");
  // -----------------------

  // all the stations without the signal should be switched from candidate --> silent
  for (auto& station : boost::make_iterator_range(sEvent.StationsBegin(), sEvent.StationsEnd())) {
    if (!station.HasRecData() || !station.GetRecData().GetTotalSignal())
      station.SetSilent();
  }

  return eSuccess;
}


VModule::ResultFlag
SdSimpleSim::Finish()
{
  if (fStatTree) {
    TFile* const file = fStatTree->GetCurrentFile();
    file->cd();
    fStatTree->Write();
    file->Close();
    delete file;
  }

  return eSuccess;
}


// ldf function from ldffinderOG
/*inline
double
SdSimpleSim::LDFFunction(const double r, const double beta, const double gamma)
{
  if (LDFFinder::fgLDFType == LDFFinder::ePL) {
    const double nearCore = kNearRadius / k1000m;
    const double relR = r / k1000m;
    return pow(relR, beta + gamma*log(relR > nearCore ? relR : nearCore));
  } else
    return pow(r/k1000m, beta)*pow((700*meter + r)/(1700*meter), beta+gamma);
}*/


double
SdSimpleSim::NKG(double N, double Rmoliere, double R, double s)
{
  if (s >= 2.25) {//Gamma(x) is only defined for x>0
    if (fVerbosity > 1) {
      ostringstream inf;
      inf <<"trying to use the NKG function outside valid range (age<=2.25): age= "<<s<<" --> using constant age = 2.24";
      WARNING(inf);
    }
    s = 2.24;
  }
  // normalize to Rmoliere
  R /= Rmoliere;
  const double C = TMath::Gamma(4.5 - s) / (TMath::Gamma(s) * TMath::Gamma(4.5 - 2*s));
  const double rho = C * N/(2*kPi*Rmoliere*Rmoliere) * pow(R, s - 2) * pow(1 + R, s - 4.5);
  return rho;
}


bool
SdSimpleSim::Noise(sevt::SEvent& sEvent,
                   const utl::TimeStamp& T0,
                   int tankID,
                   double timeOffset, double timeWindow)
{
  if (!IsNoise(timeWindow))
    return false;

  if (fVerbosity > 2)
    INFO(" ++++++++++++++++++++== generate noise!! +++++++++++++++");

  if (!sEvent.HasStation(tankID))
    sEvent.MakeStation(tankID);

  sevt::Station& tank = sEvent.GetStation(tankID);

  if (!tank.HasSimData())
    tank.MakeSimData();
  sevt::StationSimData& tankSim = tank.GetSimData();
  tankSim.SetSimulatorSignature("SdSimpleSimKG");

  /*
  // get reference coordinate system of detector (usually PampaAmarilla)
  const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
  Point pos(0,0,0,referenceCS);
  Vector dir(0,0,0,referenceCS);
  Particle fake((int)Particle::eMuon, Particle::eBackground,
                 pos, dir, TimeInterval(0), 0., 0.);
  tankSim.AddParticle(fake);
  */
  if (!tank.HasRecData())
    tank.MakeRecData();

  //sevt::StationRecData& tankRec = tank.GetRecData();
  double Signal = GenerateNoiseStation(tank, T0, timeOffset, timeWindow);

  // noise trigger type according to TOT rate (see astro-ph/0510320)

  const double timeOverThresholdRate = 1.6*hertz;

  double r = RandFlat::shoot(&fRandomEngine.GetEngine(), 0.0, 1.0);

  sevt::StationTriggerData::Algorithm tankTriggerAlgo;
  if (r > timeOverThresholdRate/fNoiseRatePerStation)
    tankTriggerAlgo = sevt::StationTriggerData::eT1Threshold;
  else
    tankTriggerAlgo = sevt::StationTriggerData::eTimeOverThreshold;

  // don't forget the tank trigger data
  if (!tank.HasTriggerData())
    tank.MakeTriggerData();

  sevt::StationTriggerData& stationTrigger = tank.GetTriggerData();

  //#warning How to properly calculate FD-T3 time offsets for the NOISE-stations?

  //int TimeOffset = 0;
  //double TimeWindow = GetCDASTriggerTimeWindow()
  stationTrigger.SetOffsetMicroSecond(timeOffset / microsecond);
  stationTrigger.SetWindowMicroSecond(timeWindow / microsecond);
  //stationTrigger.SetErrorCode(const int errorCode);
  //stationTrigger.SetPLDTrigger(const int trigger);
  stationTrigger.SetAlgorithm(tankTriggerAlgo);

  if (fVerbosity > 2) {
    ostringstream dbg;
    dbg << " noise stations: " << tankID << " signal: " << Signal << "VEM "
        << " trigger: " << (tankTriggerAlgo == sevt::StationTriggerData::eT1Threshold ? "threshold" : "TOT");
    INFO(dbg);
  }

  return true;
}


bool
SdSimpleSim::IsNoise(double TimeWindow)
{
  double NoiseEvents = 2.*TimeWindow*fNoiseRatePerStation;
  double r = RandFlat::shoot(&fRandomEngine.GetEngine(), 0.0, 1.0);

  // we support only up to one noise hit per TimeWindow .....

  if (r < NoiseEvents) {
    return true;
  }

  return false;
}


/*
  parameterized signal distribution for noise tanks, as seen from brass hybrid
  SD files (without shower stations)
*/
double
SdSimpleSim::GenerateNoiseStation(sevt::Station& tank,
                                  const utl::TimeStamp& T0,
                                  double TimeOffset,
                                  double TimeWindow)
{
  sevt::StationRecData& noiseStation = tank.GetRecData();

  double Time = RandFlat::shoot(&fRandomEngine.GetEngine(),
                                -TimeWindow, TimeWindow);

  utl::TimeStamp stationTime = T0 + utl::TimeInterval(TimeOffset+Time);

  //cout << " trig: " << T0 << endl;
  //cout << " " << Time << " " << TimeOffset << " " << TimeWindow << endl;
  //cout << "new: " << stationTime << endl;

  // parameters from analysis of data file: sd_2004_10_16_01h51.root
  // joining at mean+2*sigma !!!! = .6039
  static double Norm = 1.217e3;
  static double Mean = 4.087e-1;
  static double Sigma = 9.7614e-2;
  static double Const = 8.95;
  static double Slope = 5.37;
  static double JoiningPoint = Mean+2*Sigma;
  static double Integral_expo = exp(Const - Slope*JoiningPoint) / Slope;
  static double Integral_gaus = Norm * sqrt(kPi*2) * Sigma * 0.975;
  static double Integral = Integral_expo + Integral_gaus;
  static double p_gaus = (Integral - Integral_expo) / Integral;

  double log10VEM;
  double r0 = RandFlat::shoot(&fRandomEngine.GetEngine(), 0., 1.);
  //cout << " p_gaus: " << p_gaus << " " << r0 << endl;
  if (r0 < p_gaus) {

    // sample from gaus!
    //cout << "gauss ";
    do {
      log10VEM = RandGauss::shoot(&fRandomEngine.GetEngine(),
                                  Mean, Sigma);
    } while (log10VEM > JoiningPoint);

  } else {

    //cout << "expo ";
    // sample from expo!
    double r = RandFlat::shoot(&fRandomEngine.GetEngine(), 0., 1.);
    double N = 1./(Slope) * exp(-Slope*JoiningPoint);
    log10VEM = log(r*N*(-Slope) + exp(-Slope*JoiningPoint)) / (-Slope);
  }

  double Signal = pow(10., log10VEM);
  double OldSignal = noiseStation.GetTotalSignal();
  Signal += OldSignal;
  noiseStation.SetTotalSignal(Signal);

  const utl::TimeStamp& OldTime = noiseStation.GetSignalStartTime();

  if (OldTime>stationTime ||
      OldSignal == 0) {
    noiseStation.SetSignalStartTime(stationTime);
  }

  tank.SetHighGainSaturation(false);
  tank.SetLowGainSaturation(false);
  tank.SetCandidate();

  return Signal;
}


double
SdSimpleSim::TankIntersection(double r,
                              double phi,
                              double z,
                              double theta,
                              double tankRadius)
{
  double cosTheta = cos(theta);
  double sinTheta = sin(theta);

  double UntilBottom = 0;
  if (cosTheta != 0)
    UntilBottom = z / cos(theta);

  if (sinTheta == 0)
    return UntilBottom;

  double UntilSide = 0;
  if (sinTheta != 0)
    UntilSide = (tankRadius + r * cos(phi) -
                 (tankRadius -
                  sqrt(tankRadius*tankRadius -
                       r * r * sin(phi) * sin(phi) ))) / sin(theta);

  if (cosTheta == 0)
    return UntilSide;

  return min(UntilSide, UntilBottom);
}


double
SdSimpleSim::LTP(const double r,
                 const double /*theta*/,
                 const double /*energy*/)
{
  const double Rltp = 1;
  const double Wltp = 1;
  return 1 / (1 + exp((r - Rltp) / Wltp));
}


/*
  Electron energy distribution as parameterized by F. Nerling
 */
double
SdSimpleSim::SampleEnergy(double Emin, double Emax, double Age)
{
  Emin /= MeV;
  Emax /= MeV;

  //double a0 = 1.0;
  const double a1 = 6.42522 - 1.53183 * Age;
  const double a2 = 168.168 - 42.1368 * Age;

  double F2max = NerlingF2(Emin, a2, Age);

  double E;
  double accept;
  double r1, r2;

  // ---------INFO--------------
  int N = 0;
  // ----------------------

  do {

    // draw random number from Nerling F1
    r1 = RandFlat::shoot(&fRandomEngine.GetEngine(), 0., 1.);
    r2 = RandFlat::shoot(&fRandomEngine.GetEngine(), 0., 1.);

    E = exp(r1 * log (Emax+a1) + (1-r1) * log(Emin+a1)) - a1;

    // check with Nerling F2
    accept = NerlingF2(E, a2, Age);

    // ---------INFO--------------
    ++N;
    // ----------------------

  } while (accept/F2max < r2);

  E *= MeV;

  // ---------INFO--------------
  if (fVerbosity > 10) {
    ostringstream dbg;
    dbg << " sample-E: "
        << " Emin " << Emin
        << " Emax " << Emax
        << " age " << Age
        << E/MeV << " " << Age << " " << N;
    INFO(dbg);
  }
  // -----------------------

  return E;
}


/*
  The nerling dN/dE = F1*F2
*/
double
SdSimpleSim::NerlingF1(double E, double a1)
{
  return 1. / (E + a1);
}


double
SdSimpleSim::NerlingF2(const double E,
                       const double a2,
                       const double s)
{
  return 1 / pow((E + a2), s);
}


TabulatedFunction
SdSimpleSim::GetLogZdist(const TabulatedFunction& prof,
                         const Point& /*core*/,
                         const Vector& /*axis*/)
{
  // convert X (slant depth) to z / m
  // taking care of earth curvature
  det::Detector& Det = det::Detector::GetInstance();

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

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

  const ProfileResult& distanceVsSlantDepth =
    Det.GetAtmosphere().EvaluateDistanceVsSlantDepth();

  double atmDepthMin = distanceVsSlantDepth.MinX();
  double atmDepthMax = distanceVsSlantDepth.MaxX();

  vector <double> ValueZ;
  vector <double> ValuedNdZ;

  // one bin in depth:
  double dXslant = abs(prof[1].X() - prof[0].X());

  for (unsigned int iBin = 0; iBin < prof.GetNPoints(); ++iBin) {

    double nDmu = prof [iBin].Y();
    double xSlant = prof [iBin].X();

    double xSlantLow  = xSlant-dXslant/2;
    double xSlantHigh = xSlant+dXslant/2;

    //cout << " slant " << xSlant/g*cm*cm;

    if (xSlantLow < atmDepthMin || xSlantHigh > atmDepthMax) {
      /*
        cout << " continue "
        << " " << dXslant/g*cm*cm
        << " atmMin " << atmDepthMin/g*cm*cm
        << " atmMax " << atmDepthMax/g*cm*cm
        << endl;
      */
      continue;
    }

    const double distanceLow = distanceVsSlantDepth.Y(xSlantHigh);
    const double distanceHigh = distanceVsSlantDepth.Y(xSlantLow);
    const double distance = distanceVsSlantDepth.Y(xSlant);

    const double z = -distance;
    const double zLow = -distanceLow;
    const double zHigh = -distanceHigh;

    //cout << " profile: height " << height << " ndmu " << nDmu << flush;

    double dz = abs(zHigh - zLow);

    //cout << " dz " <<  dz << flush;

    // dont't go below ground leel
    if (z < 0) {
      //cout << endl;
      continue;
    }

    // only look at shower not before
    if (nDmu == 0) {
      //cout << endl;
      continue;
    }

    // normalisation is not required !!! ????
    double logz = log10(z);
    double dLogz = log10(dz);
    double dMuLogz = (nDmu*dXslant/dLogz) / g*cm*cm * m;

    //cout << " -> log10(z) " << logz  << " " << dMuLogz << endl;

    ValueZ.insert(ValueZ.begin(), logz);
    ValuedNdZ.insert(ValuedNdZ.begin(), dMuLogz);

  }

  /*
    for (int i=0; i< ValueZ.size(); i++ ) {

    cout << ValueZ [i] << " " << ValuedNdZ [i] << endl;
    }
  */

  /*
    TabulatedFunction result(ValueZ, ValuedNdZ);

    cout << " ---------- result -------- " << endl;
    for (int i=0; i<result.GetNPoints(); i++) {

    cout << " ->> " << result[i].X() << " " << result[i].Y() << endl;
    }
    cout << " ---------- result -------- " << endl;
    return result;
  */

  return TabulatedFunction(ValueZ, ValuedNdZ);
}


TabulatedFunction
SdSimpleSim::CalculateLogZdist(const utl::TabulatedFunction& muonProfile,
                               const double cosTheta, const double xobs)
{
  // convert X (slant depth) to z / m
  // taking care of earth curvature
  det::Detector& detector = det::Detector::GetInstance();

  const atm::ProfileResult heightProfile =
    detector.GetAtmosphere().EvaluateHeightVsDepth();
  const atm::ProfileResult depthProfile =
    detector.GetAtmosphere().EvaluateDepthVsHeight();

  const double heightXobs = heightProfile.Y(xobs);

  vector<double> valueZ;
  vector<double> valuedNdZ;

  const double atmDepthMin = depthProfile.Y(depthProfile.MaxX());
  const double atmDepthMax = depthProfile.Y(depthProfile.MinX());

  bool first = true;
  double z_old = 0;
  //  double heightX_old = 0; // unused
  double nmu_old = 0;
  for (unsigned int iBin = 0; iBin < muonProfile.GetNPoints(); ++iBin) {

    const double nmu = muonProfile[iBin].Y();
    const double xslant = muonProfile[iBin].X();
    const double xvert = xslant * cosTheta;

    if (xvert < atmDepthMin || xvert > atmDepthMax)
      continue;

    const double heightX = heightProfile.Y(xvert);
    const double z = (heightX - heightXobs) / cosTheta;

    // only look at shower not before
    if (!nmu)
      continue;

    // first bin is for init only
    if (first) {

      first = false;

      z_old = z;
      //      heightX_old = heightX; // unused
      nmu_old = nmu;
      continue;
    }

    const double dz = z - z_old;
    const double dN = nmu_old - nmu;
    const double zMean = 0.5*(z_old + z);

    // only look above the ground
    if (zMean <= 0)
      continue;

    valueZ.insert(valueZ.begin(), zMean);
    valuedNdZ.insert(valuedNdZ.begin(), dN/dz);

    z_old = z;
    //  heightX_old = heightX; // unused
    nmu_old = nmu;
  }

  return TabulatedFunction(valueZ, valuedNdZ);
}


SdSimpleSim::FitParam
SdSimpleSim::FitDecay(const TabulatedFunction& profile, bool dEdX)
{
  int nPoints = profile.GetNPoints();

  if (dEdX)
    --nPoints;

  // search for Xmax
  // start from "back" to find "last" maximum for camel-back events
  int maximumBin = 0;
  double maximum = 0;
  bool first = true;

  // use boxcar, to smooth fluctuations
  const int nBoxCar = 5;
  for (int i = 0; i < nPoints-nBoxCar; ++i) {

    // use boxcar, to smooth fluctuations
    double value = 0;
    //double depth = 0;// unused. LN.
    for (int iBoxCar = 0; iBoxCar < nBoxCar; ++iBoxCar) {

      int bin = nPoints - 1 - i - iBoxCar;
      value += profile[bin].Y();
      //depth += profile[bin].X();// unused. LN.
    }
    value /= nBoxCar;
    //depth /= nBoxCar;// unused. LN.

    if (first) {
      first = false;
      maximum = value;
      maximumBin = nPoints - 1 - i;
    } else {
      if (maximum < value) {
        maximum = value;
        maximumBin = nPoints - 1 - i;
      }
    }
  }

  // define intervall to be fit:
  // don't start at Xmax, because the exponential decay start later
  int startFitBin = maximumBin + int((nPoints - maximumBin) * 0.75);

  // ---------INFO--------------
  if (fVerbosity > 1) {
    ostringstream dbg;
    dbg << " profile fitting maxBin " << maximumBin
        << " startFitBin " << startFitBin
        << " nbins "<<nPoints;
    INFO(dbg);
  }
  // -----------------------

  if (nPoints - startFitBin < 3) {
    if (fVerbosity > 3) {
      ostringstream dbg;
      dbg << " failed profile extrapolation!\n"
             " tried to fit from bin " << startFitBin << " to " << nPoints;
      INFO(dbg);
    }

    return FitParam(0, 0);
  }

  // prepare linear-fit
  double S1 = 0;
  double Sx = 0;
  double Sxx = 0;
  double Sy = 0;
  double Sxy = 0;
  //double Syy = 0;// unused. LN.
  for (int i = startFitBin; i < nPoints; ++i) {

    const double value = log(profile[i].Y() / (dEdX ? GeV/g*cm2 : 1));
    const double depth = profile[i].X() / g*cm2;
    const double wi = 1/(.1*value * .1*value);

    S1  += wi;
    Sx  += wi * depth;
    Sxx += wi * depth * depth;
    Sy  += wi * value;
    Sxy += wi * depth * value;
    //Syy += wi * value * value;// unused. LN.
  }

  const double det = S1*Sxx - Sx*Sx;

  const double a1 = (Sxx*Sy - Sx*Sxy) / det;
  const double a2 = (-Sx*Sy + S1*Sxy) / det;

  return FitParam (a1, a2);
}


SdSimpleSim::FitParam
SdSimpleSim::FitAtmosphere(const ProfileResult& profile)
{
  // fit 10 points of the last 20th of the profile
  const double minX = profile.MinX();
  const double maxX = profile.MaxX()-50*cm; // security boundary

  const int nPoints = 10;
  const double lastPart = (maxX-minX) / 30;
  const double dX = lastPart / (nPoints-1);
  const double startX = maxX-lastPart;

  // prepare linear-fit
  double S1  = 0;
  double Sx  = 0;
  double Sxx = 0;
  double Sy  = 0;
  double Sxy = 0;
  //double Syy = 0;// unused. LN.
  for (int i = 0; i < nPoints; ++i) {

    const double x = startX + dX * i;
    const double value = log(profile.Y(x)/g*cm2);
    const double wi = 1./(.1*value*.1*value);

    S1  += wi;
    Sx  += wi * x;
    Sxx += wi * x * x;
    Sy  += wi * value;
    Sxy += wi * x * value;
    //Syy += wi * value * value;// unused. LN.

  }

  double D = S1*Sxx - Sx*Sx;

  double a1 = (Sxx*Sy - Sx * Sxy)/D;
  double a2 = (-Sx*Sy + S1*Sxy)/D;

  return FitParam (a1, a2);
}


double
SdSimpleSim::CalculateTankSignal(double coreDistance,
                                 double tankRadius, double tankHeight, double Zenith,
                                 double Age, double Rm,
                                 double Ne, double Ng, double Nmu,
                                 bool fluctuation)
{
  double cosTheta = cos (Zenith);
  double sinTheta = sin (Zenith);

  double tankProjTop = kPi * tankRadius * tankRadius * cosTheta;
  double tankProjSide =  tankHeight * 2.* tankRadius * sinTheta;
  double tankProjArea = tankProjTop + tankProjSide;


  double RhoElectrons = NKG(Ne, Rm, coreDistance, Age/fEMAgeFactor);
  double RhoGammas = NKG(Ng, Rm, coreDistance, Age/fEMAgeFactor);
  double RhoMuons = NKG(Nmu, fMuonRScale, coreDistance, Age);

  double nElectrons = RhoElectrons * tankProjArea;
  double nGammas = RhoGammas * tankProjArea;
  double nMuons = RhoMuons * tankProjArea;

  if (fVerbosity > 5) {
    ostringstream dbg;
    dbg << " +++++++ calculating tank signal at age " << Age << " and zenith " << Zenith*180./3.1415 << "deg +++++++++++\n"
           " A " << tankProjArea/m/m << "\n"
           "  coreDistance: " << coreDistance/m << "  E + G [cD/rm]: " << coreDistance/Rm << "   mu [cD/scale]: " << coreDistance/fMuonRScale << "\n"
           "   Ne : " << Ne << "   Ng: " << Ng << "   Nmu: " << Nmu << "\n"
           " rho   " << RhoElectrons << " " << RhoGammas << " " << RhoMuons << "\n"
           " nPart (no fluc) : " << nElectrons << " " << nGammas << " " << nMuons;
    INFO(dbg);
  }

  if (fluctuation) {
    nElectrons = int(RandPoisson::shoot(&fRandomEngine.GetEngine(), nElectrons));
    nGammas = int(RandPoisson::shoot(&fRandomEngine.GetEngine(), nGammas));
    nMuons = int(RandPoisson::shoot(&fRandomEngine.GetEngine(), nMuons));
  }

  fNMuons = (int)nMuons;
  fNGammas = (int)nGammas;
  fNElectrons = (int)nElectrons;

  if (fVerbosity > 5) {
    ostringstream dbg;
    dbg << " nPart " << nElectrons << " " << nGammas << " " << nMuons << endl;
    INFO(dbg);
  }

  double nParticles = nElectrons + nGammas + nMuons;

  if (nParticles <= 1) {     // no particles hit the tank
    return 0;
  }

  fSMuons = 0;
  int nSample = (fluctuation ? std::min(fNSample, int(nMuons)) : fNSample*10);
  if (nSample) {
    double Weight = nMuons / double(nSample);
    for (int iSample = 0; iSample < nSample; ++iSample) {

      const double p = RandFlat::shoot(&fRandomEngine.GetEngine(), 0, 1);

      double r;
      double phi;
      double z;
      if (p < tankProjTop/tankProjArea) { // hit tank-top

        // sample on top
        r = sqrt(RandFlat::shoot(&fRandomEngine.GetEngine(), 0, tankRadius*tankRadius));

        phi = RandFlat::shoot(&fRandomEngine.GetEngine(), 0, 2*kPi);

        z = tankHeight;

      } else {                           // hit tank-side

        // sample on side
        r = tankRadius;

        phi = asin(RandFlat::shoot(&fRandomEngine.GetEngine(), -1, 1));

        z = RandFlat::shoot(&fRandomEngine.GetEngine(), 0, tankHeight);

      }

      const double TrackLength = TankIntersection(r, phi, z, Zenith, tankRadius);

      if (fVerbosity > 10) {
        ostringstream dbg;
        dbg << " muon track/m: " << TrackLength/m;
        INFO(dbg);            
      }

      fSMuons += Weight * TrackLength / kMeterPerVEM;
    }
  } // if nMuon>0

  fSElectrons = 0;
  nSample = (fluctuation ? std::min(fNSample, int(nElectrons)) : fNSample*10);
  if (nSample) {
    double Weight = nElectrons/(double)nSample;
    for (int iSample = 0; iSample < nSample; ++iSample) {

      //double Eem = MeanEnergy;
      static double Emin = 1e-2 * MeV;
      static double Emax = 1e+5 * MeV;
      double Eem = SampleEnergy(Emin, Emax, Age);

      fSElectrons += Weight * Eem * kTrackPerGEV / kMeterPerVEM;
    }
  } // if nElectrons>0

  fSGammas = 0;
  nSample = (fluctuation ? std::min(fNSample, int(nGammas)) : fNSample*10);
  if (nSample) {
    double Weight = nGammas/(double)nSample;
    for (int iSample = 0; iSample < nSample; ++iSample) {

      //double Eem = MeanEnergy;
      static double Emin = 1e-2 * MeV;
      static double Emax = 1e+5 * MeV;
      double Eem = SampleEnergy(Emin, Emax, Age);

      fSGammas += Weight * Eem  * kTrackPerGEV / kMeterPerVEM;
    }
  } // if nGammas>0

  return fSMuons + fSElectrons + fSGammas;
}


/*
  See GAP note 2005-059
*/

static double fgTriggerThetaBins [5] = { 15, 25, 35, 45, 55 };
static double fgTriggerS1000Bins [4] = { 4.5, 13.5, 40.5, 121.5 };

static double fgTriggerPropShalf [5] [4] = {
  { 1.9, 2.2, 3.6, 5.1 },
  { 1.9, 3.0, 3.5, 4.8 },
  { 2.5, 3.3, 4.1, 5.1 },
  { 3.5, 3.8, 4.2, 5.4 },
  { 4.2, 4.8, 4.8, 5.0 }
};

static double fgTriggerPropShalfError [5] [4] = {
  { .6, .9, .4, .2 },
  { .9, .5, .5, .3 },
  { .7, .6, .5, .2 },
  { .4, .5, .5, .2 },
  { .9, .4, .4, .2 }
};

double
SdSimpleSim::T1TriggerProbability(double signal, double S1000, double theta)
{
  static double n = 2.8; // or 3

  if (S1000 > fgTriggerS1000Bins[3]) {
    S1000 = fgTriggerS1000Bins[3];
  }

  if (S1000 < fgTriggerS1000Bins[0]) {
    S1000 = fgTriggerS1000Bins[0];
  }

  if (theta > fgTriggerThetaBins[4]) {
    theta = fgTriggerThetaBins[4];
  }

  if (theta < fgTriggerThetaBins[0]) {
    theta = fgTriggerThetaBins[0];
  }

  // ---------INFO--------------
  if (fVerbosity > 3) {
    ostringstream dbg;
    dbg << " pT1: " << theta << " " << S1000 << " " << signal;
    INFO(dbg);
  }
  // -----------------------

  int thetaBin;
  for (thetaBin = 0; thetaBin < 4; ++thetaBin) {
    if (theta >= fgTriggerThetaBins[thetaBin] && theta <= fgTriggerThetaBins[thetaBin+1]) {
      break;
    }
  }

  int s1000Bin;
  for (s1000Bin = 0; s1000Bin < 3; ++s1000Bin) {
    if (S1000 >= fgTriggerS1000Bins[s1000Bin] && S1000 <= fgTriggerS1000Bins[s1000Bin+1]) {
      break;
    }
  }

  double S_half_11 = fgTriggerPropShalf[thetaBin][s1000Bin];
  double S_half_12 = fgTriggerPropShalf[thetaBin][s1000Bin+1];
  double S_half_21 = fgTriggerPropShalf[thetaBin+1][s1000Bin];
  double S_half_22 = fgTriggerPropShalf[thetaBin+1][s1000Bin+1];

  double S_half_error_11 = fgTriggerPropShalfError[thetaBin][s1000Bin];
  double S_half_error_12 = fgTriggerPropShalfError[thetaBin][s1000Bin+1];
  double S_half_error_21 = fgTriggerPropShalfError[thetaBin+1][s1000Bin];
  double S_half_error_22 = fgTriggerPropShalfError[thetaBin+1][s1000Bin+1];

  // extrapolate in theta direction

  double delta_s1000 = (S1000 - fgTriggerS1000Bins[s1000Bin]) / (fgTriggerS1000Bins[s1000Bin+1] - fgTriggerS1000Bins[s1000Bin]);
  double S_half_1 = S_half_11 + (S_half_12 - S_half_11) * delta_s1000;
  double S_half_2 = S_half_21 + (S_half_22 - S_half_21) * delta_s1000;
  double S_half_1_error = S_half_error_11 + (S_half_error_12 - S_half_error_11) * delta_s1000;
  double S_half_2_error = S_half_error_21 + (S_half_error_22 - S_half_error_21) * delta_s1000;

  double delta_theta = (theta - fgTriggerThetaBins[thetaBin]) / (fgTriggerThetaBins[thetaBin+1] - fgTriggerThetaBins[thetaBin]);
  double S_half = S_half_1 + (S_half_2 - S_half_1) * delta_theta;
  double S_half_error = S_half_1_error + (S_half_2_error - S_half_1_error) * delta_theta;

  // ---------INFO--------------
  if (fVerbosity > 3) {
    ostringstream dbg;
    dbg << " S_0.5 " << S_half << " " << S_half_error;
    INFO(dbg);
  }
  // -----------------------

  double p = pow(signal, n);
  p /= (p + pow(S_half, n));
  return p;
}


double
SdSimpleSim::GetCDASTriggerTimeWindow(const TimeStamp& time)
{
  const double switchTriggerGPS = 795.3e6; // GPS sec, as taken from plotting data vs. time

  if (time.GetGPSSecond() < switchTriggerGPS)
    return 120. * microsecond;
  else
    return 20. * microsecond;
}