G4TankSimulatorOG/G4TankSimulator.cc

Go to the documentation of this file.

#include "G4TankSimulator.h"
#include "G4TankPhysicsListCustomization.h"
#include "G4TankPhysicsList.h"
#include "G4TankConstruction.h"
#include "G4TankPrimaryGenerator.h"
#include "G4TankTrackingAction.h"
#include "G4TankStackingAction.h"
#include "G4TankSteppingAction.h"
#include "G4TankEventAction.h"
#include "G4TankRunAction.h"
#include "G4TankVisManager.h"
#include "G4TankFastCerenkov.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 <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>

#include <CLHEP/Random/Random.h>

#include <tls/Geant4Manager.h>

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


const sdet::Station* G4TankSimulator::fCurrentDetectorStation = 0;
SEvent::StationIterator G4TankSimulator::fCurrentEventStationIt;
StationSimData::ParticleIterator G4TankSimulator::fCurrentParticleIt;

bool G4TankSimulator::fgMuCapture;


G4TankSimulator::G4TankSimulator() :
  fRunManager(0),
  fUImanager(0),
  fVisManager(0),
  fGeoVisOn(0),
  fTrajVisOn(0),
  fUseGlobalPhysicsList(false),
  fSignalSeparationMode(eStandard),
  fDetectorConstructed(0),
  fFastMode(false),
  fEventId("")
{ }


G4TankSimulator::~G4TankSimulator()
{ }


VModule::ResultFlag
G4TankSimulator::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("G4TankSimulator");
  Branch visB = topB.GetChild("visualization");

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

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

  Branch muCaptureB = topB.GetChild("muCapture");
  muCaptureB.GetData(fgMuCapture);

  Branch globalPhysicsListB = topB.GetChild("globalPhysicsList");
  globalPhysicsListB.GetData(fUseGlobalPhysicsList);

  string sepMode;
  topB.GetChild("signalSeparationMode").GetData(sepMode);
  if (sepMode == "Standard")
    fSignalSeparationMode = eStandard;
  else if (sepMode == "Universality")
    fSignalSeparationMode = eUniversality;

  if ( (fGeoVisOn || fTrajVisOn) && !fVisManager) {
    fVisManager = new G4TankVisManager;
  }
  if (!fRunManager) {
    if (fUseGlobalPhysicsList)
      fRunManager = tls::Geant4Manager::GetInstance().GetRunManager();
    else
      fRunManager = new G4RunManager;
  }

  fUImanager = G4UImanager::GetUIpointer();

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

  if (fUseGlobalPhysicsList) {
    INFO("using global PhysicsList from Geant4Manager");
    tls::Geant4Manager::GetInstance().AddVisManager(fVisManager);

    G4TankPhysicsListCustomization* physicsListCustomization =
      new G4TankPhysicsListCustomization(fFastMode);

    INFO("Constructing G4Tank");
    const SDetector& sDetector = Detector::GetInstance().GetSDetector();
    const SDetector::StationIterator sIt = sDetector.AllStationsBegin();
    fCurrentDetectorStation = &(*sIt);
    fgTankConstruction = new G4TankConstruction(); //reads fCurrentDetectorStation

    tls::Geant4Customization custom("G4TankSimulator", fgTankConstruction, physicsListCustomization);
    custom.SetPrimaryGenerator(new G4TankPrimaryGenerator);
    fStackingAction = new G4TankStackingAction();
    custom.SetStackingAction(fStackingAction);
    custom.SetTrackingAction(new G4TankTrackingAction);
    custom.SetSteppingAction(new G4TankSteppingAction);
    custom.SetEventAction(new G4TankEventAction);
    custom.SetRunAction(new G4TankRunAction);

    tls::Geant4Manager::GetInstance().AddCustomization(custom);
  }
  else {
    INFO("using PhysicsList from G4TankPhysicsList");
    fRunManager->SetUserInitialization(new G4TankPhysicsList(fFastMode));

    fRunManager->SetUserAction(new G4TankPrimaryGenerator);

    fStackingAction = new G4TankStackingAction();
    fRunManager->SetUserAction(fStackingAction);
    fRunManager->SetUserAction(new G4TankTrackingAction);
    fRunManager->SetUserAction(new G4TankSteppingAction);
    fRunManager->SetUserAction(new G4TankEventAction);
    fRunManager->SetUserAction(new G4TankRunAction);

    fUImanager->ApplyCommand("/run/verbose 0");
    fUImanager->ApplyCommand("/event/verbose 0");
    fUImanager->ApplyCommand("/tracking/verbose 0");
  }

  fDetectorConstructed = false;

  return eSuccess;
}


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

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

  SEvent& sEvent = event.GetSEvent();

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

    if (!fUseGlobalPhysicsList) {
      INFO("Constructing G4Tank");

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

      const SDetector& sDetector = Detector::GetInstance().GetSDetector();
      fCurrentDetectorStation = &sDetector.GetStation(sIt->GetId());

      // Construct detector geometry and register w/ RunManager
      fRunManager->SetUserInitialization(new G4TankConstruction);

      fRunManager->Initialize();
    }

    fDetectorConstructed = true;

    if (fGeoVisOn || fTrajVisOn) {
       fVisManager->Initialize();
       fUImanager->ApplyCommand("/vis/open VRML2FILE");
       fUImanager->ApplyCommand("/vis/scene/create");
       fUImanager->ApplyCommand("/vis/sceneHandler/attach");
       fUImanager->ApplyCommand("/vis/scene/add/volume");
       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");
    }

  }

  if (fUseGlobalPhysicsList)
    tls::Geant4Manager::GetInstance().Customize("G4TankSimulator");

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


VModule::ResultFlag
G4TankSimulator::RunFull(Event& event)
{
  // Get the SEvent
  SEvent& sEvent = event.GetSEvent();

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

  for (SEvent::StationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt) {

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

    if (sIt->HasSimData()) {

      StationSimData& simData = sIt->GetSimData();

      // full-blown tank sim (no optimization of photon tracking)
      simData.SetSimulatorSignature("G4TankSimulatorFullOG");

      const int stationId = sIt->GetId();

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

      fCurrentEventStationIt = sIt;
      fCurrentDetectorStation = &sDetector.GetStation(stationId);
      fStackingAction->SetStation(fCurrentDetectorStation);

      ConstructTraces(*sIt);

      for (StationSimData::ParticleIterator pIt = simData.ParticlesBegin();
           pIt != simData.ParticlesEnd(); ++pIt) {

        fCurrentParticleIt = pIt;

        fRunManager->BeamOn(1);

      }
    }
  }

  return eSuccess;
}


VModule::ResultFlag
G4TankSimulator::RunFast(Event& event)
{
  SEvent& sEvent = event.GetSEvent();
  const SDetector& sDetector = Detector::GetInstance().GetSDetector();

  // find first station with SimData
  for (fCurrentEventStationIt = sEvent.StationsBegin();
       fCurrentEventStationIt != sEvent.StationsEnd() &&
       !fCurrentEventStationIt->HasSimData(); ++fCurrentEventStationIt) ;

  if (fCurrentEventStationIt == sEvent.StationsEnd()) {
    INFO("No stations with SimData found.");
    return eContinueLoop;
  }

  const int id = fCurrentEventStationIt->GetId();
  fCurrentDetectorStation = &sDetector.GetStation(id);
  if (!fDetectorConstructed) {
    fgTankConstruction = new G4TankConstruction();
    fRunManager->SetUserInitialization(fgTankConstruction);
    fRunManager->Initialize();
    fDetectorConstructed = true;
  }

  for (SEvent::StationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt) {

    if (!sIt->HasSimData())
      continue;

    StationSimData& simData = sIt->GetSimData();

    simData.SetSimulatorSignature("G4TankSimulatorOG");

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

    const int stationId = sIt->GetId();

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

    fCurrentEventStationIt = sIt;
    fCurrentDetectorStation = &sDetector.GetStation(stationId);

    ConstructTraces(*sIt);

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

    for (StationSimData::ParticleIterator pIt = simData.ParticlesBegin();
         pIt != simData.ParticlesEnd(); ++pIt) {

      fCurrentParticleIt = pIt;

      fRunManager->BeamOn(1);

    } // particles

  } // stations

  G4cerr << endl;

  return eSuccess;
}


VModule::ResultFlag
G4TankSimulator::Finish()
{
  if (fUseGlobalPhysicsList)
    tls::Geant4Manager::GetInstance().NotifyDelete();
  else {
    delete fRunManager;
    fRunManager = 0;
  }
  return eSuccess;
}


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

  // Add pmt sim data only if the station contains particles
  if (sSim.ParticlesEnd() == sSim.ParticlesBegin())
    return;

  for (sevt::Station::PMTIterator pmtIt = station.PMTsBegin();
       pmtIt != station.PMTsEnd(); ++pmtIt)
    if (!pmtIt->HasSimData())
      pmtIt->MakeSimData();
}


void
G4TankSimulator::AddPhoton(const int nPMT, const double peTime)
  const
{
  PMTSimData& pmtSim = fCurrentEventStationIt->GetPMT(nPMT).GetSimData();

  const int type = fCurrentParticleIt->GetType();
  const int source = fCurrentParticleIt->GetSource();

  sevt::StationConstants::SignalComponent component = sevt::StationConstants::eTotal;
  sevt::StationConstants::SignalComponent 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);
  }
}