CachedDirectInjector.cc

Go to the documentation of this file.

#include <utl/config.h>

#include <cmath>
#include <sstream>
#include <vector>
#include <fstream>

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <io/CorsikaUtilities.h>

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

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

#include <det/Detector.h>

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

#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/GeometryUtilities.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/TimeStamp.h>
#include <utl/Math.h>
#include <utl/Reader.h>
#include <utl/TabularStream.h>
#include <utl/NucleusProperties.h>

#include "CachedDirectInjector.h"

using namespace CachedDirectInjectorOG;
using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace utl;
using namespace std;


namespace CachedDirectInjectorOG {

  istream&
  operator>>(istream& is, InjectedParticle& p)
  {
    string ppp;
    return
      is >> ws >> ppp
         >> ws >> p.fId
         >> ws >> p.fR
         >> ws >> p.fPhi
         >> ws >> p.fHeight
         >> ws >> p.fSTime
         >> ws >> p.fType
         >> ws >> p.fSource
         >> ws >> p.fPosX >> ws >> p.fPosY >> ws >> p.fPosZ
         >> ws >> p.fDirX >> ws >> p.fDirY >> ws >> p.fDirZ
         >> ws >> p.fTime
         >> ws >> p.fWeight
         >> ws >> p.fEk
         >> ws;
  }

  inline
  double
  PlaneFrontTime(const CoordinateSystemPtr showerCS, const Point& corePosition,
                 const Point& position)
  {
    return (corePosition.GetZ(showerCS) - position.GetZ(showerCS)) / kSpeedOfLight;
  }

  inline
  double
  Round(const double div, const double val)
  {
    return round(div*val)/div;
  }


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

    AttributeMap useAtt;
    useAtt["use"] = string("yes");

    if (topB) {

      //Check for file with additional particles to be injected
      Branch partTrackB = topB.GetChild("ParticleInjectionFile", useAtt);
      if (partTrackB) {
        const string filename = partTrackB.Get<string>();
        ifstream file(filename);
        if (!file.is_open()) {
          ERROR("Could not open particle-injection file.");
          return eFailure;
        } else {
          ostringstream info;
          info << "Successfully opened particle injection file " << filename;
          INFO(info);

          string line;
          while (getline(file, line)) {
            const InjectedParticle particle = boost::lexical_cast<InjectedParticle>(line);
            fTankToInjectedParticles[particle.fId].push_back(particle);
          }
        }
      }

      //Check for particle-tank distance and other option
      Branch particleTankDistB = topB.GetChild("ParticleTankDistanceCut");
      if (particleTankDistB)
        particleTankDistB.GetData(fParticleTankDistanceCut);
      Branch maxPart = topB.GetChild("MaxNumberOfParticles");
      if (maxPart)
        maxPart.GetData(fMaxParticles);

      //Check for dense stations option
      Branch denseStations = topB.GetChild("UseDenseStations");
      if (denseStations)
        denseStations.GetData(fUseDense);
      if(fUseDense){
        Branch denseStationsOnly = topB.GetChild("UseDenseStationsOnly");
        if (denseStationsOnly)
          denseStationsOnly.GetData(fUseDenseOnly);
      }

    }
    else{
      ERROR("Cannot find XML scheme for direct injection.");
      return eFailure;
    }

    return eSuccess;
  }


  VModule::ResultFlag
  CachedDirectInjector::Run(evt::Event& event)
  {

    ostringstream info;

    if (!event.HasSimShower()) {
      ERROR("Current event does not have a simulated shower.");
      return eFailure;
    }

    const ShowerSimData& simShower = event.GetSimShower();
    const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
    const Point& corePosition = simShower.GetPosition();
    const TimeStamp& coreTime = simShower.GetTimeStamp();

    if (!event.HasSEvent())
      event.MakeSEvent();
    SEvent& sEvent = event.GetSEvent();

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

    vector<const sdet::Station*> dStationVec;
    vector<sevt::Station*> stationVec;

    vector<double> sRShowerFrameVec;
    vector<double> sPhiShowerFrameVec;

    for (sdet::SDetector::StationIterator sdIt = sDetector.StationsBegin();
         sdIt != sDetector.StationsEnd(); ++sdIt) {

      const sdet::Station& dStation = *sdIt;
      const int sId = dStation.GetId();
      const CoordinateSystemPtr dStationCS = dStation.GetLocalCoordinateSystem();

      if (!fUseDense && dStation.IsDense())
        continue; // only NOT dense stations

      if (fUseDenseOnly && !dStation.IsDense())
        continue; // only dense stations

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

      Station& station = sEvent.GetStation(sId);

      if (!station.HasSimData())
        station.MakeSimData();
      StationSimData& sSim = station.GetSimData();

      if (station.HasScintillator() && !station.GetScintillator().HasSimData())
        station.GetScintillator().MakeSimData();

      sSim.ClearParticleList();

      const Point sPos = sdIt->GetPosition();
      const double sR = sPos.GetRho(showerCS);

      const TimeStamp planeTime =
        coreTime + TimeInterval(PlaneFrontTime(showerCS, corePosition, sPos));
      sSim.SetPlaneFrontTime(planeTime);

      dStationVec.push_back(&dStation);
      stationVec.push_back(&station);
      sRShowerFrameVec.push_back(sR);
      sPhiShowerFrameVec.push_back(sPos.GetPhi(showerCS));

      // let's inject some dummies, so we can, for example, force a trigger
      // in the not dense stantions ...
      if (!dStation.IsDense() && fTankToInjectedParticles.size()) {
        if (fTankToInjectedParticles.count(sId)) {
          for (const InjectedParticle& p : fTankToInjectedParticles[sId]) {
            const Vector myDir(p.fDirX, p.fDirY, p.fDirZ, dStationCS);
            const Point pShiftedPosition(p.fR, p.fPhi, p.fHeight, dStationCS, Point::kCylindrical);
            const Particle myParticle(p.fType, (utl::Particle::Source)p.fSource,
                                      pShiftedPosition, myDir, p.fSTime, p.fWeight, p.fEk);
            station.AddParticle(myParticle);
          }
        }
      }
    }

    const unsigned int nStations = dStationVec.size();
    fEarliestTime.clear();

    unsigned int nParticles = 0;
    unsigned int nParticlesInjected = 0;

    for (fParticleIt = simShower.GroundParticlesBegin();
         fParticleIt != simShower.GroundParticlesEnd(); ++fParticleIt) {

      //if (nParticles % 1000 == 0)
      //std::cout << "Particles : " << nParticles << '\r';

      const Point pPosition = fParticleIt->GetPosition();
      const double pRShowerFrame = pPosition.GetRho(showerCS);
      const double pPhiShowerFrame = pPosition.GetPhi(showerCS);
      const double pTime = fParticleIt->GetTime().GetInterval();

      // do not inject unknown particles
      const int particleType = io::Corsika::CorsikaToPDG(fParticleIt->GetType());
      if (particleType == utl::Particle::eUndefined)
        continue;

      // Check particle time relative to the shower front
      // (see CachedShowerGenerator.cc for explanation)
      const double pPlaneFrontDelay =
        pTime - PlaneFrontTime(showerCS, corePosition, pPosition);
        if (pPlaneFrontDelay < -0.1*ns)
          continue;

      for (unsigned int si = 0; si < nStations; ++si) {

        const sdet::Station& dStation = *dStationVec[si];
        Station& station = *stationVec[si];
        const unsigned int sId = dStation.GetId();
        const CoordinateSystemPtr dStationCS = dStation.GetLocalCoordinateSystem();

        // Check if particle is nearer than chosen distance
        const double dPhi = NormalizeAngleMinusPiPi(pPhiShowerFrame - sPhiShowerFrameVec[si]);
        if (std::abs(sRShowerFrameVec[si] - pRShowerFrame) > fParticleTankDistanceCut ||
            std::abs(dPhi)*sRShowerFrameVec[si] > fParticleTankDistanceCut)
          continue;

        TimeInterval& earliestTime = fEarliestTime[sId];
        if (earliestTime.GetInterval() < numeric_limits<double>::min())
          earliestTime = TimeInterval::Max();

        // Is particle already inside station ?
        if (dStation.IsInsideStation(fParticleIt->GetPosition()))
           continue;

        bool particleHit = false;

        if (dStation.HasScintillator())
          particleHit = dStation.GetScintillator().IsHit(fParticleIt->GetPosition(),
                                                         fParticleIt->GetDirection());
        if (!particleHit)
          particleHit = dStation.IsHit(fParticleIt->GetPosition(),
                                       fParticleIt->GetDirection());

        if (particleHit) {

          Particle newParticle(*fParticleIt);

          const double sThickness = dStation.GetThickness();
          const double sHeight = dStation.GetHeight() + 2*sThickness;
          double injectionHeight = sHeight;
          if (dStation.HasScintillator()) {
            const double scintHeight = dStation.GetScintillator().GetMaxHeight();
            injectionHeight = max(injectionHeight, scintHeight);
          }
          const Line particleLine(fParticleIt->GetPosition(), fParticleIt->GetDirection());
          const Plane topPlane(Point(0, 0, injectionHeight*utl::meter, dStationCS),
                               Vector(0, 0, 1, dStationCS));
          Point hitPoint = Intersection(topPlane, particleLine);
          const double timeShift =
            pTime + (fParticleIt->GetPosition() - hitPoint).GetMag() / kSpeedOfLight;

          newParticle.SetPosition(hitPoint);
          newParticle.SetTime(timeShift);

          if (timeShift < earliestTime)
            earliestTime = timeShift;

          station.AddParticle(newParticle);
          ++nParticlesInjected;

          info.str("");
          info << "\nInjecting particle with ";
          info << "energy = " << newParticle.GetKineticEnergy()/utl::GeV << " GeV ";
          info << "at time = "   << newParticle.GetTime()/utl::ns << " ns \n";
          info << "in (x,y,z) = (" << newParticle.GetPosition().GetX(dStationCS)/utl::meter << ", "
               << newParticle.GetPosition().GetY(dStationCS)/utl::meter << ", "
               << newParticle.GetPosition().GetZ(dStationCS)/utl::meter << ") m\n";
          //INFO(info);

        } // if particle hits a detector

      } // end cycle over stations

      if (nParticlesInjected >= fMaxParticles) {
        info.str("");
        info << "Reached max number of " << nParticlesInjected << " injected particles.";
        INFO(info);
        return eSuccess;
      }

      ++nParticles;
    }

    info.str("");
    info << nParticlesInjected << " particles injected over "
         << nParticles << " particles processed.";
    INFO(info);

    for (unsigned si = 0; si < nStations; ++si) {

      const unsigned sId = dStationVec[si]->GetId();
      const TimeStamp newTime = coreTime + fEarliestTime[sId];
      StationSimData& sSim = stationVec[si]->GetSimData();
      const TimeStamp sTime = sSim.GetPlaneFrontTime();

      if (!sTime || newTime < sTime)
        sSim.SetPlaneFrontTime(newTime);

    }

    return eSuccess;
  }


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

}