G4StationSimulator.cc

Go to the documentation of this file.

#include "G4StationSimulator.h"
#include "G4StationPhysicsList.h"
#include "G4StationConstruction.h"
#include "G4StationPrimaryGenerator.h"
#include "G4StationTrackingAction.h"
#include "G4StationStackingAction.h"
#include "G4StationSteppingAction.h"
#include "G4StationEventAction.h"
#include "G4StationRunAction.h"
#include "G4StationVisManager.h"
#include "G4StationFastCerenkov.h"
#include "G4OutputHandler.h"

#include <G4RunManager.hh>
#include <G4UImanager.hh>
#include <G4VisManager.hh>

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

#include <det/Detector.h>

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

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

#include <mdet/MDetector.h>
#include <mdet/Counter.h>

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

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/Particle.h>
#include <utl/ParticleCases.h>
#include <utl/ShowerParticleIterator.h>
#include <utl/TimeDistribution.h>
#include <utl/TimeDistributionAlgorithm.h>

#include <cstddef>
#include <iostream>
#include <sstream>

//using namespace utl;
using namespace fwk;
using namespace std;
using namespace det;
using namespace evt;
using namespace sdet;
using namespace sevt;
using namespace mdet;
using namespace mevt;
using namespace cevt;


namespace G4StationSimulatorOG {

  bool G4StationSimulator::Setup::fgLock = false;

  G4StationSimulator::Current G4StationSimulator::fgCurrent;
  bool G4StationSimulator::fgMuCapture = false;


  VModule::ResultFlag
  G4StationSimulator::Init()
  {
    auto& cc = *CentralConfig::GetInstance();
    const utl::Branch topB = cc.GetTopBranch("G4StationSimulator");
    const utl::Branch visB = topB.GetChild("visualization");

    visB.GetChild("geometry").GetData(fGeoVisOn);
    visB.GetChild("trajectories").GetData(fTrajVisOn);
    visB.GetChild("renderfile").GetData(fRenderFile);

    const utl::Branch verB = topB.GetChild("verbosity");
    if (verB)
      verB.GetData(fVerbosity);

    const utl::Branch trackModeB = topB.GetChild("fullTrackMode");
    trackModeB.GetData(fTrackMode);

    const utl::Branch fastModeB = topB.GetChild("fastMode");
    fastModeB.GetData(fFastMode);

    const utl::Branch rangeCutsB = topB.GetChild("rangeCuts");
    if (rangeCutsB) {
      rangeCutsB.GetChild("default").GetData(fRangeCutDefault);
      rangeCutsB.GetChild("muonplus").GetData(fRangeCutMuonPlus);
      rangeCutsB.GetChild("muonminus").GetData(fRangeCutMuonMinus);
      rangeCutsB.GetChild("electron").GetData(fRangeCutElectron);
      rangeCutsB.GetChild("positron").GetData(fRangeCutPositron);
      rangeCutsB.GetChild("gamma").GetData(fRangeCutGamma);
      rangeCutsB.GetChild("opticalphoton").GetData(fRangeCutOpticalPhoton);
    }

    const utl::Branch martaB = cc.GetTopBranch("CachedShowerRegenerator").GetChild("MARTA");
    // Needed to be consistent with CachedShowerRegenerator area of injection (NB: variable loaded with AugerUnits)
    // MARTA needs to be checked before AMIGA because of its smaller simulation radius
    if (martaB) {
      martaB.GetChild("Simulate").GetData(fMARTAEnabled);
      if (fMARTAEnabled) {
        INFO("MARTA RPC G4 simulation has been activated");
        martaB.GetChild("Radius").GetData(fSimulationRadius);
        // Auger to CLHEP unit conversion
        fSimulationRadius *= CLHEP::m / utl::m;
      }
    }

    const utl::Branch umdB = cc.GetTopBranch("CachedShowerRegenerator").GetChild("UMD");
    // Needed to be consistent with CachedShowerRegenerator area of injection (NB: variable loaded with AugerUnits)
    if (umdB) {
      umdB.GetChild("Simulate").GetData(fUMDEnabled);
      if (fUMDEnabled) {
        INFO("UMD tank+ground G4 simulation has been activated");
        umdB.GetChild("MaxRadius").GetData(fSimulationRadius);
        fSimulationRadius += 1*utl::cm;
        // Auger to CLHEP unit conversion
        fSimulationRadius *= CLHEP::m / utl::m;
        // Parameters for underground propagation
        const utl::Branch ugrdB = topB.GetChild("underGroundSim");
        ugrdB.GetChild("umd_scinti_yield").GetData(fScintYield);
        fScintYield /= utl::keV;  // unit-less number
        // Tracking or not the optical photons in umd scintillators
        ugrdB.GetChild("umdFastMode").GetData(fUMDFastMode);
        // Tracking or not low energetic paricles underground
        ugrdB.GetChild("umdFastProp").GetData(fUMDFastProp);
        ugrdB.GetChild("onlyMuons").GetData(fUMDOnlyMuons);
        const utl::Branch cutsB = ugrdB.GetChild("energyLimits");
        cutsB.GetChild("muons").GetData(fUMDMuonsCut);
        cutsB.GetChild("electrons").GetData(fUMDElectCut);
        cutsB.GetChild("gammas").GetData(fUMDGammaCut);
        cutsB.GetChild("other").GetData(fUMDOtherCut);
        fUMDMuonsCut *=  CLHEP::eV / utl::eV;
        fUMDElectCut *=  CLHEP::eV / utl::eV;
        fUMDGammaCut *=  CLHEP::eV / utl::eV;
        fUMDOtherCut *=  CLHEP::eV / utl::eV;
      }
    }
    const utl::Branch muCaptureB = topB.GetChild("muCapture");
    muCaptureB.GetData(fgMuCapture);

    const auto sepMode = topB.GetChild("signalSeparationMode").Get<string>();
    if (sepMode == "Standard")
      fSignalSeparationMode = eStandard;
    else if (sepMode == "Universality")
      fSignalSeparationMode = eUniversality;

    if ((fGeoVisOn || fTrajVisOn) && !fVisManager)
      fVisManager = new G4StationVisManager;

    if (!fRunManager)
      fRunManager = new G4RunManager;

    fUImanager = G4UImanager::GetUIpointer();

    G4OutputHandler* const logger = new G4OutputHandler();
    fUImanager->SetCoutDestination(logger);

    fRunManager->SetUserInitialization(new G4StationPhysicsList(*this, fFastMode));

    fRunManager->SetUserAction(new G4StationPrimaryGenerator);

    fStackingAction = new G4StationStackingAction();
    fRunManager->SetUserAction(fStackingAction);
    fRunManager->SetUserAction(new G4StationTrackingAction(fUMDEnabled));
    fRunManager->SetUserAction(new G4StationSteppingAction(fMARTAEnabled));
    fRunManager->SetUserAction(new G4StationEventAction(fUMDEnabled));
    fRunManager->SetUserAction(new G4StationRunAction);

    switch (fVerbosity) {
    case 1:
      fUImanager->ApplyCommand("/run/verbose 1");
      fUImanager->ApplyCommand("/event/verbose 0");
      fUImanager->ApplyCommand("/tracking/verbose 0");
      break;
    case 2:
      fUImanager->ApplyCommand("/run/verbose 1");
      fUImanager->ApplyCommand("/event/verbose 1");
      fUImanager->ApplyCommand("/tracking/verbose 0");
      break;
    case 3:
      fUImanager->ApplyCommand("/run/verbose 1");
      fUImanager->ApplyCommand("/event/verbose 1");
      fUImanager->ApplyCommand("/tracking/verbose 1");
      break;
    default:
      fUImanager->ApplyCommand("/run/verbose 0");
      fUImanager->ApplyCommand("/event/verbose 0");
      fUImanager->ApplyCommand("/tracking/verbose 0");
    }

    fDetectorConstructed = false;

    return eSuccess;
  }


  VModule::ResultFlag
  G4StationSimulator::Run(Event& event)
  {
    if (!event.HasSEvent()) {
      ERROR("SEvent does not exist.");
      return eFailure;
    }

    if (fMARTAEnabled && !event.HasCEvent()) {
      event.MakeCEvent();
      //WARNING("CEvent does not exist. Making it now. This should not be a problem if you are using the ParticleInjector module.");
    }

    G4Random::setTheEngine(&RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine());

    SEvent& sEvent = event.GetSEvent();

    if (!fDetectorConstructed && sEvent.GetNumberOfStations() > 0) {

      INFO("Constructing G4Station");

      const auto sIt = sEvent.StationsBegin();
      if (sIt == sEvent.StationsEnd()) {
        ERROR("No stations!");
        return eFailure;
      }

      Setup setup(Detector::GetInstance().GetSDetector().GetStation(sIt->GetId()));

      // Construct detector geometry and register w/ RunManager
      fRunManager->SetUserInitialization(new G4StationConstruction(fSimulationRadius, fUMDEnabled, fScintYield, fMARTAEnabled));

      fRunManager->Initialize();

      fDetectorConstructed = true;

      if (fGeoVisOn || fTrajVisOn) {
        fVisManager->Initialize();
        fUImanager->ApplyCommand(("/vis/open " + fRenderFile).c_str());
        fUImanager->ApplyCommand("/vis/scene/create");
        fUImanager->ApplyCommand("/vis/sceneHandler/attach");
        fUImanager->ApplyCommand("/vis/scene/add/volume");
        fUImanager->ApplyCommand("/vis/viewer/set/viewpointThetaPhi 0. 0.");
        fUImanager->ApplyCommand("/vis/viewer/set/targetPoint 0 0 0");
        fUImanager->ApplyCommand("/vis/viewer/zoom 1");
        fUImanager->ApplyCommand("/vis/viewero/set/style/wireframe");
        fUImanager->ApplyCommand("/vis/drawVolume");
        fUImanager->ApplyCommand("/vis/scene/notifyHandlers");
        fUImanager->ApplyCommand("/vis/viewer/update");
      }
      if (fTrajVisOn) {
        fUImanager->ApplyCommand("/tracking/storeTrajectory 1");
        fUImanager->ApplyCommand("/vis/scene/add/trajectories");
      }

    }

    return fFastMode ? RunFast(event) : RunFull(event);
  }


  VModule::ResultFlag
  G4StationSimulator::RunFull(Event& event)
  {
    // full-blown tank sim (no optimization of photon tracking) with
    // G4-internal tracking of all particles (including photons)
    if (fTrackMode)
      fSimulatorSignature = "G4StationSimulatorFullTrackOG";
    else {
      // full-blown tank sim (no optimization of photon tracking) but with
      // photon killing according to the PMT quantum efficiency
      fSimulatorSignature = "G4StationSimulatorFullOG";
    }

    INFO(fSimulatorSignature);

    // Get the SEvent
    SEvent& sEvent = event.GetSEvent();

    if (!event.HasMEvent() && fUMDEnabled)
      event.MakeMEvent();

    // Get the SDetector
    const SDetector& sDetector = Detector::GetInstance().GetSDetector();

    for (auto& station : sEvent.StationsRange()) {

      if (!fDetectorConstructed) {
        ERROR("The detector hasn't been initialized when looping over stations");
        return eFailure;
      }

      const int stationId = station.GetId();

      // Add Stations and SimData to CEvent (MARTA)
      if (fMARTAEnabled) {
        CEvent& cEvent = event.GetCEvent();
        if (!cEvent.HasStation(stationId)) {
          cEvent.MakeStation(stationId);
          cEvent.GetStation(stationId).MakeSimData();
        }
      }

      if (station.HasSimData()) {
        sevt::StationSimData& simData = station.GetSimData();
        simData.SetSimulatorSignature(fSimulatorSignature);

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

        if (numParticles)
          G4cerr << stationId << ':' << numParticles << ' ' << flush;
        else
          continue;

        /* 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);

        Setup setup;
        setup.Set(station);
        setup.Set(sDetector.GetStation(stationId));

        if (fUMDEnabled && fgCurrent.GetDetectorStation().ExistsAssociatedCounter()) {
          const auto counterId = fgCurrent.GetDetectorStation().GetAssociatedCounterId();

          MEvent& mEvent = event.GetMEvent();
          if (!mEvent.HasCounter(counterId))
            mEvent.MakeCounter(counterId);
          setup.Set(mEvent.GetCounter(counterId));

          const MDetector& mDetector = Detector::GetInstance().GetMDetector();
          setup.Set(mDetector.GetCounter(counterId));
        }

        ConstructTraces(station);

        for (auto& particle : simData.ParticlesRange()) {
          setup.Set(particle);
          fRunManager->BeamOn(1);
          // Add particles from simulation to CEvent data (particles in the RPC's)
          if (fMARTAEnabled) {
            CEvent& cEvent = event.GetCEvent();
            cevt::StationSimData& martaSimData = cEvent.GetStation(stationId).GetSimData();
            FillRPCSimData(martaSimData);
          }
        }
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  G4StationSimulator::RunFast(Event& event)
  {
    fSimulatorSignature = "G4StationSimulatorOG";
    INFO(fSimulatorSignature);

    // Get the SEvent
    SEvent& sEvent = event.GetSEvent();

    if (!event.HasMEvent() && fUMDEnabled)
      event.MakeMEvent();

    // Get the SDetector
    const SDetector& sDetector = Detector::GetInstance().GetSDetector();

    if (!fDetectorConstructed) {
      // find first station with SimData
      sevt::Station* station = nullptr;
      for (auto& s : sEvent.StationsRange())
        if (s.HasSimData()) {
          station = &s;
          break;
        }
      if (!station) {
        INFO("No stations with SimData found.");
        return eSuccess;
      }

      Setup setup(sDetector.GetStation(*station));
      fStationConstruction = new G4StationConstruction(fSimulationRadius, fUMDEnabled, fScintYield, fMARTAEnabled);
      fRunManager->SetUserInitialization(fStationConstruction);
      fRunManager->Initialize();
      fDetectorConstructed = true;
    }

    for (auto& station : sEvent.StationsRange()) {

      const int stationId = station.GetId();

      // Add Stations and SimData to CEvent (MARTA)
      if (fMARTAEnabled) {
        CEvent& cEvent = event.GetCEvent();
        if (!cEvent.HasStation(stationId)) {
          cEvent.MakeStation(stationId);
          cEvent.GetStation(stationId).MakeSimData();
        }
      }

      if (!station.HasSimData())
        continue;

      sevt::StationSimData& simData = station.GetSimData();

      simData.SetSimulatorSignature(fSimulatorSignature);

      const unsigned long numParticles = simData.GetNParticles();

      if (numParticles)
        cerr << "StationId " << stationId << " : Particles " << numParticles << endl;
      else
        continue;

      /* 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);

      Setup setup;
      setup.Set(station);
      setup.Set(sDetector.GetStation(station));

      if (fUMDEnabled) {
        if (fgCurrent.GetDetectorStation().ExistsAssociatedCounter()) {
          const auto counterId = fgCurrent.GetDetectorStation().GetAssociatedCounterId();
          MEvent& mEvent = event.GetMEvent();
          if (!mEvent.HasCounter(counterId))
            mEvent.MakeCounter(counterId);
          setup.Set(mEvent.GetCounter(counterId));
          const MDetector& mDetector = Detector::GetInstance().GetMDetector();
          setup.Set(mDetector.GetCounter(counterId));
        }
      }

      ConstructTraces(station);

      // prod the Cerenkov process update it's tank description data
      G4StationFastCerenkov::GetDataFromConstruction();

      for (auto& particle : simData.ParticlesRange()) {
        setup.Set(particle);
        fRunManager->BeamOn(1);
        // Add particles from simulation to CEvent data
        if (fMARTAEnabled) {
          CEvent& cEvent = event.GetCEvent();
          cevt::StationSimData& martaSimData = cEvent.GetStation(stationId).GetSimData();
          FillRPCSimData(martaSimData);
        }
      }

    }

    // flush all remaining G4 output (otherwise its flushed after the program ends)
    G4cout << endl;

    return eSuccess;
  }


  VModule::ResultFlag
  G4StationSimulator::Finish()
  {
    delete fRunManager;
    fRunManager = nullptr;

    return eSuccess;
  }


  void
  G4StationSimulator::ConstructTraces(sevt::Station& station)
    const
  {
    const sevt::StationSimData& sSim = station.GetSimData();

    // Add pmt sim data only if the station contains particles
    if (!sSim.GetNParticles())
      return;

    for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType))
      if (!pmt.HasSimData())
        pmt.MakeSimData();
  }


  void
  G4StationSimulator::AddInjectedParticle(const size_t modId, const size_t stripId, const utl::Particle& particle)
    const
  {
    if (!fgCurrent.GetEventUMDCounter().HasSimData())
      fgCurrent.GetEventUMDCounter().MakeSimData();

    if (!fgCurrent.GetEventUMDCounter().HasModule(modId))
      fgCurrent.GetEventUMDCounter().MakeModule(modId);

    mevt::Module& module = fgCurrent.GetEventUMDCounter().GetModule(modId);

    if (!module.HasScintillator(stripId))
      module.MakeScintillator(stripId);

    mevt::Scintillator& scinti = module.GetScintillator(stripId);

    if (!scinti.HasSimData())
      scinti.MakeSimData();

    mevt::ScintillatorSimData& simData = scinti.GetSimData();

    // Load particle in the proper event scintillator
    // (NOTE: any particle entering in more than one strip (or bouncing back to same strip) is stored in each corresponding
    // scintillator strip but with different type:
    // 1) as utl::Particle::eShower for the first strip (or first time in a strip) and
    // 2) as utl::Particle::eCornerClipping in the following stirp(s)
    simData.AddParticle(particle);

    const auto type = particle.GetType();
    if (type == utl::Particle::eMuon || type == utl::Particle::eAntiMuon) {
      if (particle.GetSource() != utl::Particle::eCornerClipping)
        simData.IncrementNumberOfInjectedMuons();
      else
        simData.IncrementNumberOfCornerClippingMuons();
    } else if (type == utl::Particle::eElectron || type == utl::Particle::ePositron)
      simData.IncrementNumberOfElectrons();
  }


  void
  G4StationSimulator::AddUMDPhoton(const size_t modId, const size_t stripId, const double peTime)
    const
  {
    // peTime is in Auger units!

    auto& counter = fgCurrent.GetEventUMDCounter();
    if (!counter.HasSimData())
      counter.MakeSimData();

    if (!counter.HasModule(modId))
      counter.MakeModule(modId);

    mevt::Module& module = counter.GetModule(modId);

    if (!module.HasScintillator(stripId))
      module.MakeScintillator(stripId);

    mevt::Scintillator& scinti = module.GetScintillator(stripId);

    if (!scinti.HasSimData())
      scinti.MakeSimData();

    mevt::ScintillatorSimData& simData = scinti.GetSimData();

    simData.AddPhotonTime(peTime);  // Auger units
  }


  void
  G4StationSimulator::AddPhoton(const int nPMT, const double peTime)
    const
  {
    /* In principle, we should print an ERROR and throw a non existent component
      exception. Once the hardware creation/destruction/update issues have been
      resolved within geant, the warning below can be changed to an error and an
      exception thrown. */
    if (!fgCurrent.GetEventStation().HasPMT(nPMT)) {
      WARNING("Attempt to add photoelectron to PMT that does not exist. Doing nothing.");
      return;
    }

    PMTSimData& pmtSim = fgCurrent.GetEventStation().GetPMT(nPMT).GetSimData();

    const auto type = fgCurrent.GetParticle().GetType();
    const auto source = fgCurrent.GetParticle().GetSource();

    auto component = sevt::StationConstants::eTotal;
    auto extraComponent = sevt::StationConstants::eTotal;

    switch (fSignalSeparationMode) {
    case eStandard:
      switch (type) {
      case OFFLINE_ELECTRONS:
        component = sevt::StationConstants::eElectron;
        if (source == utl::Particle::eShowerFromMuonDecay)
          extraComponent = sevt::StationConstants::eShowerMuonDecayElectron;
        break;
      case OFFLINE_PHOTON:
        component = sevt::StationConstants::ePhoton;
        if (source == utl::Particle::eShowerFromMuonDecay)
          extraComponent = sevt::StationConstants::eShowerMuonDecayPhoton;
        break;
      case OFFLINE_MUONS:
        component = sevt::StationConstants::eMuon;
        break;
      case OFFLINE_HADRONS:
        component = sevt::StationConstants::eHadron;
        break;
      default:
        cerr << "Unknown particle type" << endl;
        exit(-1);
      }
      break;
    case eUniversality:
      switch (type) {
      case OFFLINE_ELECTRONS:
        if (source == utl::Particle::eShowerFromLocalHadron)
          component = sevt::StationConstants::eShowerLocalHadronElectron;
        else if (source == utl::Particle::eShowerFromMuonDecay)
          component = sevt::StationConstants::eShowerMuonDecayElectron;
        else
          component = sevt::StationConstants::eElectron;
        break;
      case OFFLINE_PHOTON:
        if (source == utl::Particle::eShowerFromLocalHadron)
          component = sevt::StationConstants::eShowerLocalHadronPhoton;
        else if (source == utl::Particle::eShowerFromMuonDecay)
          component = sevt::StationConstants::eShowerMuonDecayPhoton;
        else
          component = sevt::StationConstants::ePhoton;
        break;
      case OFFLINE_MUONS:
        component = sevt::StationConstants::eMuon;
        break;
      case OFFLINE_HADRONS:
        component = sevt::StationConstants::eHadron;
        break;
      default:
        cerr << "Unknown particle type" << endl;
        exit(-1);
      }
      break;
    }

    if (!pmtSim.HasPETimeDistribution())
      pmtSim.MakePETimeDistribution();
    pmtSim.GetPETimeDistribution().AddTime(peTime);

    if (!pmtSim.HasPETimeDistribution(component))
      pmtSim.MakePETimeDistribution(component);
    pmtSim.GetPETimeDistribution(component).AddTime(peTime);

    if (extraComponent != sevt::StationConstants::eTotal) {
      if (!pmtSim.HasPETimeDistribution(extraComponent))
        pmtSim.MakePETimeDistribution(extraComponent);
      pmtSim.GetPETimeDistribution(extraComponent).AddTime(peTime);
    }
  }


  void
  G4StationSimulator::FillRPCSimData(cevt::StationSimData& simData)
    const
  {
    auto& rpcParticleListMap = G4StationSteppingAction::fgRPCParticleListMap;
    for (const auto& rpcParticleList : rpcParticleListMap) {
      const auto& particleList = rpcParticleList.second;
      for (const auto& p : particleList)
        simData.AddParticle(p);
    }
    rpcParticleListMap.clear();

    G4StationSteppingAction::fgRPCTrackIDListMap.clear();
  }

}