FastTankSimulator.cc

Go to the documentation of this file.

#include "FastTankSimulator.h"

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

#include <det/Detector.h>

#include <evt/Event.h>

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

#include <sevt/PMT.h>
#include <sevt/PMTSimData.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>

#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/RandomEngine.h>
#include <utl/TimeDistribution.h>
#include <utl/Math.h>
#include <utl/Reader.h>

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

using namespace FastTankSimulatorOG;
using namespace fwk;
using namespace det;
using namespace evt;
using namespace sdet;
using namespace sevt;
using namespace utl;
using namespace std;


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


FastTankSimulator::FastTankSimulator() :
  fInitialParticle(
    utl::Particle::eUndefined,
    utl::Particle::eShower,
    utl::Point(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
    utl::Vector(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
    utl::TimeInterval(0),
    0
  )
{ }


VModule::ResultFlag
FastTankSimulator::Init()
{
  Branch topBranch =
    CentralConfig::GetInstance()->GetTopBranch("FastTankSimulator");

  topBranch.GetChild("dEdXMuon").GetData(fdEdXMuon);
  topBranch.GetChild("dEdXElectron").GetData(fdEdXElectron);
  topBranch.GetChild("CerenkovMinMuon").GetData(fCerenkovMinMuon);
  topBranch.GetChild("CerenkovMinElectron").GetData(fCerenkovMinElectron);
  topBranch.GetChild("EMinPhoton").GetData(fEMinPhoton);
  topBranch.GetChild("PhotoElectronRate").GetData(fPhotoElectronRate);
  topBranch.GetChild("WaterRIndex").GetData(fWaterRIndex);
  topBranch.GetChild("PECumulativeProb").GetData(fPECumulativeProbability);

  vector<double> photonInteractionLength;
  topBranch.GetChild("PhotonInteractionLength").GetData(photonInteractionLength);
  vector<double> pairProductionProbability;
  topBranch.GetChild("PairProductionProbability").GetData(pairProductionProbability);
  vector<double> energyTable;
  topBranch.GetChild("EnergyTable").GetData(energyTable);

  fPairProductionProbability =
    utl::TabulatedFunction(energyTable, pairProductionProbability);
  fPhotonInteractionLength =
    utl::TabulatedFunction(energyTable, photonInteractionLength);

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

  return eSuccess;
}


VModule::ResultFlag
FastTankSimulator::Run(evt::Event& theEvent)
{
  if (!theEvent.HasSEvent()) {
    ERROR("Non-existent SEvent.");
    return eFailure;
  }

  sevt::SEvent& theSEvent = theEvent.GetSEvent();

  const sdet::SDetector& theSDetector =
    det::Detector::GetInstance().GetSDetector();

  for (sevt::SEvent::StationIterator sIt = theSEvent.StationsBegin();
       sIt != theSEvent.StationsEnd(); ++sIt)

    if (sIt->HasSimData()) {

      sevt::StationSimData& simData = sIt->GetSimData();
      simData.SetSimulatorSignature("FastTankSimulatorOG"); // should this necessarily be same as module name?
      fCurrentEventStation = &(*sIt);

      fCurrentDetectorStation = &theSDetector.GetStation(sIt->GetId());

      fTankHeight = fCurrentDetectorStation->GetHeight();

      const double tankRadius = fCurrentDetectorStation->GetRadius();

      fTankRadiusSquared = Sqr(tankRadius);

      const unsigned long numParticles =
        simData.ParticlesEnd() - simData.ParticlesBegin();

      /*
       * Counting the number of particles that are actually simulated (post
       * station-level thinning). Due to resampling simulation option that
       * uses cycles with maximum particle number per cycle, we always check
       * if particles have already been simulated.
       */
      const unsigned int nParticlesAlreadySimulated = simData.GetTotalSimParticleCount();
      if (!nParticlesAlreadySimulated)
        simData.SetTotalSimParticleCount(numParticles);
      else
        simData.SetTotalSimParticleCount(nParticlesAlreadySimulated + numParticles);

      for (sevt::StationSimData::ParticleIterator pIt = sIt->GetSimData().ParticlesBegin();
           pIt != sIt->GetSimData().ParticlesEnd(); ++pIt) {

        fInitialParticle = *pIt;

        switch (fInitialParticle.GetType()) {
        case utl::Particle::ePositron:
        case utl::Particle::eElectron:
          SimulateElectrons(fInitialParticle);
          break;
        case utl::Particle::ePhoton:
          SimulatePhotons(fInitialParticle);
          break;
        case utl::Particle::eMuon:
        case utl::Particle::eAntiMuon:
          SimulateMuons(fInitialParticle);
          break;
        }

      }

    }

  return eSuccess;
}


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


double
FastTankSimulator::CalculateDistanceInTank(const Particle& particle)
{
  const CoordinateSystemPtr cs =
    fCurrentDetectorStation->GetLocalCoordinateSystem();

  double x, y, z;
  boost::tie(x, y, z) = particle.GetPosition().GetCoordinates(cs);
  double pX, pY, pZ;
  boost::tie(pX, pY, pZ) = particle.GetDirection().GetCoordinates(cs);

  const double a = particle.GetDirection().GetRho(cs);
  const double b = x*pX + y*pY;
  const double b2 = Sqr(b);
  const double c = particle.GetPosition().GetRho2(cs) - fTankRadiusSquared;
  const double ac = a*c;

  // Vertical component
  const double vertDistance = (pZ < 0) ? -z/pZ : (fTankHeight - z)/pZ;

  // Horizontal component
  if (a) {
    const double horizDistance = (b2 >= ac) ? (-b + sqrt(b2 - ac))/a : 0;
    return min(vertDistance, horizDistance);
  }

  return vertDistance;
}


double
FastTankSimulator::CalculateIntegratedCerenkovRate(const Particle& particle,
                                                   const double energyLoss)
{
  double ratio;
  double mass;

  switch (particle.GetType()) {

  case utl::Particle::eElectron:
  case utl::Particle::ePositron:
    mass = utl::kElectronMass;
    ratio = fPhotoElectronRate/fdEdXElectron;
    break;

  case utl::Particle::eMuon:
  case utl::Particle::eAntiMuon:
    mass = utl::kMuonMass;
    ratio = fPhotoElectronRate/fdEdXMuon;
    break;

  default:
    return 0;

  }

  const double energy = particle.GetKineticEnergy();
  const double energyDifference = energy - energyLoss;

  if (energyDifference > 0)
    return ratio*energyDifference *
      (1 - (mass/(2*energyDifference*(Sqr(fWaterRIndex) - 1))) *
       log((energy/energyLoss) *
           (energyLoss + 2*mass)/(energy + 2*mass)));

  return 0;
}


void
FastTankSimulator::CalculatePhotoElectrons(const Particle& particle,
                                           const double cerenkovRate)
{
  for (sevt::Station::PMTIterator pmtIt = fCurrentEventStation->PMTsBegin();
       pmtIt != fCurrentEventStation->PMTsEnd(); ++pmtIt) {

    if (!pmtIt->HasSimData())
      pmtIt->MakeSimData();

    sevt::PMTSimData& thePMTSimData = pmtIt->GetSimData();

    const unsigned int numberPhotoElectrons =
      RandPoisson::shoot(&fRandomEngine->GetEngine(), cerenkovRate);

    for (unsigned int countPE = 0; countPE < numberPhotoElectrons; ++countPE) {

      // For each PE, choose a time slot according to the
      // tabulated probability

      const double randNum =
        RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);

      // Guaranteed to terminate
      int j;
      for (j = 0; randNum > fPECumulativeProbability[j]; ++j)
        ;

      const double time = particle.GetTime().GetInterval() + j*ns;

      if (!thePMTSimData.HasPETimeDistribution())
        thePMTSimData.MakePETimeDistribution();

      pmtIt->GetSimData().GetPETimeDistribution().AddTime(time);

      sevt::StationConstants::SignalComponent component = sevt::StationConstants::eTotal;

      switch (fInitialParticle.GetType()) {
      case utl::Particle::eElectron:
      case utl::Particle::ePositron:
        component = sevt::StationConstants::eElectron;
        break;
      case utl::Particle::ePhoton:
        component = sevt::StationConstants::ePhoton;
        break;
      case utl::Particle::eMuon:
      case utl::Particle::eAntiMuon:
        component = sevt::StationConstants::eMuon;
        break;
      }

      if (component != sevt::StationConstants::eTotal) {
        if (!thePMTSimData.HasPETimeDistribution(component))
          thePMTSimData.MakePETimeDistribution(component);
        thePMTSimData.GetPETimeDistribution(component).AddTime(time);
      }

    }

  }

}


void
FastTankSimulator::SimulateElectrons(const Particle& particle)
{
  const double energy = particle.GetKineticEnergy();

  if (energy < fCerenkovMinElectron)
    return;

  const double distance = CalculateDistanceInTank(particle);

  const double energyLoss =
    max(energy - fdEdXElectron*distance, fCerenkovMinElectron);

  const double rate = CalculateIntegratedCerenkovRate(particle, energyLoss);
  CalculatePhotoElectrons(particle, rate);
}


void
FastTankSimulator::SimulateMuons(const Particle& particle)
{
  const double energy = particle.GetKineticEnergy();

  if (energy < fCerenkovMinMuon)
    return;

  const double distance = CalculateDistanceInTank(particle);

  const double energyLoss =
    max(energy - fdEdXMuon*distance, fCerenkovMinMuon);

  const double rate = CalculateIntegratedCerenkovRate(particle, energyLoss);
  CalculatePhotoElectrons(particle, rate);
}


void
FastTankSimulator::SimulatePhotons(const Particle& particle)
{
  const double energy = particle.GetKineticEnergy();

  // If energy is below 100 keV, return without any further simulation. The old
  // way of calculating the index into the table caused problems on Bruce's
  // laptop. It was a retarded, and un-safe, way of calculating the proper
  // index into the interaction length table, anyway. So, we now use a
  // pre-calculated energy table, and have some tighter sanity checks before
  // continuing with the photon simulation.

  if (energy < fPhotonInteractionLength.GetX(0))
    return;

  //if (energy < fEnergyTable.front())
    //return;

  //unsigned int index = 0;

  //for ( ; index < fEnergyTable.size(); ++index)
    //if (energy <= fEnergyTable.at(index))
      //break;

  //const unsigned int index = min(max(0, int(lookupValue)), int(fPhotonInteractionLength.size() - 1));

  // We assume that if the energy is greater than 10 GeV, the interaction
  // length, and pair production probability are constant, and set at the
  // values obtained for 10 GeV

  // Use this for both the interaction length and production probability size
  // since they are over the same size energy table

  const unsigned int nPoints = fPhotonInteractionLength.GetNPoints();

  const double interactionLength =
    (energy > fPhotonInteractionLength.GetX(nPoints-1)) ?
      fPhotonInteractionLength.GetY(nPoints-1) : fPhotonInteractionLength.Y(energy);

  const double pathLength =
    -log(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1)) *
    interactionLength;

  const double distanceInTank = CalculateDistanceInTank(particle);

  if (pathLength < distanceInTank) {

    const Point pos =
      particle.GetPosition() + pathLength*particle.GetDirection();

    const TimeInterval time =
      particle.GetTime() + utl::TimeInterval(pathLength/kSpeedOfLight);

    const Particle newParticle(
      particle.GetType(),
      particle.GetSource(),
      pos,
      particle.GetDirection(),
      time,
      particle.GetWeight(),
      energy
    );

    const double randNum = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);

    const double productionProbability =
      (energy > fPairProductionProbability.GetX(nPoints-1)) ?
        fPairProductionProbability.GetY(nPoints-1) : fPairProductionProbability.Y(energy);

    if (randNum < productionProbability)
      SimulatePairProduction(newParticle);
    else
      SimulateComptonScattering(newParticle);

  }
}


void
FastTankSimulator::SimulatePairProduction(const Particle& particle)
{
  utl::Particle newElectron(
    utl::Particle::eElectron,
    particle.GetSource(),
    particle.GetPosition(),
    particle.GetDirection(),
    particle.GetTime(),
    particle.GetWeight(),
    particle.GetKineticEnergy()
  );

  utl::Particle newPositron(
    utl::Particle::ePositron,
    particle.GetSource(),
    particle.GetPosition(),
    particle.GetDirection(),
    particle.GetTime(),
    particle.GetWeight(),
    particle.GetKineticEnergy()
  );

  const double maxKineticEnergy = particle.GetKineticEnergy() - 2*kElectronMass;
  const double kineticEnergy =
    RandFlat::shoot(&fRandomEngine->GetEngine(), 0, maxKineticEnergy);

  newElectron.SetKineticEnergy(kineticEnergy);
  newPositron.SetKineticEnergy(maxKineticEnergy - kineticEnergy);

  SimulateElectrons(newElectron);
  SimulateElectrons(newPositron);
}


void
FastTankSimulator::SimulateComptonScattering(const Particle& particle)
{
  const double energy = particle.GetKineticEnergy();

  // Minimum energy for final photon
  const double energyMin = 1/(1/energy + 2/kElectronMass);

  // Simplified KN distribution with rejection
  const double b = energyMin/energy + energy/energyMin;
  double energyFinal;
  double a;

  do {

    energyFinal =
      RandFlat::shoot(&fRandomEngine->GetEngine(), energyMin, energy);

    a = energyFinal/energy + energy/energyFinal;

  } while (RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1) > a/b);

  // Scattering angles (final photon wrt initial photon)
  const double cosTheta = 1 + kElectronMass * (1/energy - 1/energyFinal);
  const double sinTheta = sqrt(1 - Sqr(cosTheta));

  const double phi =
    RandFlat::shoot(&fRandomEngine->GetEngine(), 0, kTwoPi);
  const double cosPhi = cos(phi);
  const double sinPhi = sin(phi);

  // Construct two unit vectors perpendicular to the initial photon
  const CoordinateSystemPtr cs =
    fCurrentDetectorStation->GetLocalCoordinateSystem();

  const Vector dir = particle.GetDirection();

  double pX, pY, pZ;
  boost::tie(pX, pY, pZ) = dir.GetCoordinates(cs);

  const double uR = dir.GetRho(cs);

  double vX, vY;
  double wX, wY, wZ;

  if (uR > 0) {

    vX = -pY/uR;
    vY = pX/uR;
    wX = -pZ*vY;
    wY = pZ*vX;
    wZ = pX*vY - pY*vX;

  } else {

    vX = 1;
    vY = 0;
    wX = 0;
    wY = 1;
    wZ = 0;

  }

  const Vector u = cosTheta * dir +
    sinTheta * Vector(cosPhi*vX+sinPhi*wX, cosPhi*vY+sinPhi*wY, sinPhi*wZ, cs);

  Vector v = energy * dir - energyFinal * u;

  v.Normalize();

  // This is the final photon
  utl::Particle newPhoton(
    utl::Particle::ePhoton,
    particle.GetSource(),
    particle.GetPosition(),
    u,
    particle.GetTime(),
    particle.GetWeight(),
    energyFinal
  );

  // This is the final electron
  const double electronEnergy = energy - energyFinal;
  utl::Particle newElectron(
    utl::Particle::eElectron,
    particle.GetSource(),
    particle.GetPosition(),
    v,
    particle.GetTime(),
    particle.GetWeight(),
    electronEnergy
  );

  // Do the simulation
  SimulateElectrons(newElectron);
  SimulatePhotons(newPhoton);
}