G4XTankSimulator.cc

Go to the documentation of this file.

#include "G4XTankSimulator.h"
#include "G4XTankPhysicsListCustomization.h"
#include "G4XTankPhysicsList.h"
#include "G4XTankConstruction.h"
#include "G4XTankPrimaryGenerator.h"
#include "G4XTankTrackingAction.h"
#include "G4XTankStackingAction.h"
#include "G4XTankSteppingAction.h"
#include "G4XTankEventAction.h"
#include "G4XTankRunAction.h"
#include "G4XTankVisManager.h"
#include "G4XTankFastCerenkov.h"
#include "G4XTankPMT.h"
#include "G4XOutputHandler.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 <mevt/MEvent.h>
#include <mevt/CounterSimData.h>
#include <mdet/MDetector.h>
#include <mdet/Counter.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/Particle.h>
#include <utl/ShowerParticleIterator.h>
#include <utl/TimeDistribution.h>
#include <utl/TimeDistributionAlgorithm.h>
#include <utl/ParticleCases.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 mevt;
using namespace mdet;
using namespace evt;
using namespace det;
using namespace G4XTankSimulatorAG;


G4XTankPMT* G4XTankSimulator::fgPMTAction = nullptr;
const sdet::Station* G4XTankSimulator::fCurrentDetectorStation = nullptr;
SEvent::StationIterator G4XTankSimulator::fCurrentEventStationIt;
StationSimData::ParticleIterator G4XTankSimulator::fCurrentParticleIt;

bool G4XTankSimulator::fgMuCapture = false;


sevt::StationConstants::SignalComponent
G4XTankSimulator::GetComponentId(const StationSimData::ConstParticleIterator curParticle)
{
  switch (curParticle->GetType()) {
  case OFFLINE_ELECTRONS:
    return sevt::StationConstants::eElectron;
  case OFFLINE_PHOTON:
    return sevt::StationConstants::ePhoton;
  case OFFLINE_MUONS:
    return sevt::StationConstants::eMuon;
  case OFFLINE_HADRONS:
    return sevt::StationConstants::eHadron;
  default:
    return sevt::StationConstants::eTotal;
  }
}


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

  AttributeMap useAtt;
  useAtt["use"] = "yes";

  Branch OnGrdB = topB.GetChild("storeOnGrdParticles", useAtt);
  if (OnGrdB) {
    fStoreOnGrdPart = true;
    cout << "Particles on ground will be stored in Counter particle list structure X-XX\n";
  }

  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");
  if (globalPhysicsListB)
    globalPhysicsListB.GetData(fUseGlobalPhysicsList);

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

  if (!fRunManager) {
    if (fUseGlobalPhysicsList)
      fRunManager = tls::Geant4Manager::GetInstance().GetRunManager();
    else
      fRunManager = new G4RunManager;
  }

  fUImanager = G4UImanager::GetUIpointer();

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

  //Setting up counters' depth
  const SDetector& sDetector = Detector::GetInstance().GetSDetector();
  SDetector::StationIterator sIt = sDetector.StationsBegin();
  ostringstream ss;

  while (sIt != sDetector.StationsEnd()) {

    if (sIt->ExistsAssociatedCounter()) {
      fCurrentDetectorStation = &(*sIt);
      ss.str("");
      ss << "Station " << fCurrentDetectorStation->GetId() << " has associated counter";
      INFO(ss);
      break;
    }

    ++sIt;
  }

  if (!fCurrentDetectorStation) {
    ERROR("Could not find any station with an associated muon counter");
    return eFailure;
  }

  const MDetector& mDetector = det::Detector::GetInstance().GetMDetector();
  const mdet::Counter& currentCounter = mDetector.GetCounter(fCurrentDetectorStation->GetAssociatedCounterId());
  const mdet::Module& module = *(currentCounter.ModulesBegin());

  const CoordinateSystemPtr dStationCS = fCurrentDetectorStation->GetLocalCoordinateSystem();
  double CtrDepth = -1 * currentCounter.GetPosition().GetZ(dStationCS) / utl::m*CLHEP::m;
  double ModDepth = -1 * module.GetPosition().GetZ(dStationCS) / utl::m*CLHEP::m;

  ss.str("");
  ss << "\n"
        "Counter depth: " << CtrDepth/CLHEP::m << " m\n"
        "Module depth: " << ModDepth/CLHEP::m << " m\n"
        "-> Using " << (CtrDepth>ModDepth ? "Counter" : "Module") << " depth";
  INFO(ss);

  fDepth = CtrDepth > ModDepth ? CtrDepth : ModDepth;

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

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

    G4XTankPhysicsListCustomization* const physicsListCustomization = new G4XTankPhysicsListCustomization(fFastMode);

    INFO("Constructing G4XTank");

    fgTankConstruction = new G4XTankConstruction( fDepth ); //reads fCurrentDetectorStation

    tls::Geant4Customization custom("G4XTankSimulator", fgTankConstruction, physicsListCustomization);
    custom.SetPrimaryGenerator(new G4XTankPrimaryGenerator);
    fStackingAction = new G4XTankStackingAction();
    custom.SetStackingAction(fStackingAction);
    custom.SetTrackingAction(new G4XTankTrackingAction);
    fSteppingAction = new G4XTankSteppingAction();
    custom.SetSteppingAction(fSteppingAction);
    custom.SetEventAction(new G4XTankEventAction);
    custom.SetRunAction(new G4XTankRunAction);

    tls::Geant4Manager::GetInstance().AddCustomization(custom);

  } else {

    INFO("using PhysicsList from G4XTankPhysicsList");
    fRunManager->SetUserInitialization(new G4XTankPhysicsList(fFastMode));

    fRunManager->SetUserAction(new G4XTankPrimaryGenerator);

    fStackingAction = new G4XTankStackingAction();
    fRunManager->SetUserAction(fStackingAction);
    fRunManager->SetUserAction(new G4XTankTrackingAction);
    fSteppingAction = new G4XTankSteppingAction();
    fRunManager->SetUserAction(fSteppingAction);
    fRunManager->SetUserAction(new G4XTankEventAction);
    fRunManager->SetUserAction(new G4XTankRunAction);

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

  fDetectorConstructed = false;

  return eSuccess;
}


VModule::ResultFlag
G4XTankSimulator::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 G4XTank");

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

      fgTankConstruction = new G4XTankConstruction(fDepth); // reads fCurrentDetectorStation

      // Construct detector geometry and register w/ RunManager
      fRunManager->SetUserInitialization(fgTankConstruction); // reads fCurrentDetectorStation

      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("G4XTankSimulator");

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


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

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

  // Get MEvent
  if (!event.HasMEvent())
    event.MakeMEvent();

  MEvent& mEvent = event.GetMEvent();

  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 << "On ground -> StatId " << stationId << ':' << numParticles << " particle(s)" << flush;
      else
        continue;

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

      ConstructTraces(*sIt);

      //Check what happened in the associated muon counter
      if (!fCurrentDetectorStation->ExistsAssociatedCounter())
        continue;

      // Get the associated muon counter ID
      const int mId = fCurrentDetectorStation->GetAssociatedCounterId();

      // Get the counter
      if (!mEvent.HasCounter(mId))
        mEvent.MakeCounter(mId);

      mevt::Counter& eCounter = mEvent.GetCounter(mId);

      // Get the counterSimData
      if (!eCounter.HasSimData())
        eCounter.MakeSimData();

      CounterSimData& cSim = eCounter.GetSimData();

      // Get the muon counter depth
      fSteppingAction->fCounterDepth = -fDepth;

      // Run the particle propagation
      for (StationSimData::ParticleIterator pIt = simData.ParticlesBegin();
           pIt != simData.ParticlesEnd(); ++pIt) {

        fCurrentParticleIt = pIt;
        smartPointer parent(new Particle(*fCurrentParticleIt));
        fSteppingAction->SetImpingingParticle(parent);

        if (fStoreOnGrdPart) {
          Particle part = *fCurrentParticleIt;
          cSim.AddGrdParticle(part);
        }

        fRunManager->BeamOn(1);

      }

      cout << fSteppingAction->fPartColl.size()  << " particles underground at the counter's depth "
           << fSteppingAction->fCounterDepth/CLHEP::m << " m"
              " on station number " << sIt->GetId() << " ";

      const CoordinateSystemPtr CounterCS = fCurrentDetectorStation->GetLocalCoordinateSystem();
      G4XTankSteppingAction::PartCollectionIt pCollIt;

      unsigned int muons_ugrd = 0;
      //Fill the CounterSimData with the particles which reached the counters depth
      for (pCollIt = fSteppingAction->fPartColl.begin(); pCollIt != fSteppingAction->fPartColl.end(); ++pCollIt) {
        PartData pData = (*pCollIt);
        /*
        cout << "pData.Pname "   << pData.Pname   << " - " << pData.KinE/CLHEP::MeV << " (MeV)" << endl;
        cout << "pData.PDGCode " << pData.PDGCode << " - " << pData.Time/CLHEP::ns  << " (ns)" << endl;
        cout << "pData.PParent->GetType() " << pData.PParent->GetType() << " - " << (pData.PParent->GetType()==utl::Particle::eMuon||pData.PParent->GetType()==utl::Particle::eAntiMuon?"SHOWER MUON!!!!":"OTHER PARTICLE")<< endl;
        cout << "pData.TrackID " << pData.TrackID << endl;

        cout << "pData.x " << pData.x/CLHEP::m << endl;
        cout << "pData.y " << pData.y/CLHEP::m << endl;
        cout << "pData.z " << pData.z/CLHEP::m << endl;

        cout << "pData.px " << pData.px << endl;
        cout << "pData.py " << pData.py << endl;
        cout << "pData.pz " << pData.pz << endl;
        */
        if (pData.PParent->GetType() == utl::Particle::eMuon ||
            pData.PParent->GetType() == utl::Particle::eAntiMuon)
          ++muons_ugrd;

        Point pPos(pData.x / CLHEP::m * utl::m, pData.y / CLHEP::m * utl::m, pData.z / CLHEP::m * utl::m, CounterCS);
        Vector pDir(pData.px, pData.py, pData.pz, CounterCS);
        //Creates new particle
        Particle part(pData.PDGCode, utl::Particle::eShower, pPos, pDir, pData.Time / CLHEP::ns * utl::ns,
                      1, pData.KinE/CLHEP::MeV * utl::MeV );
        part.SetParent(*pData.PParent);
        cSim.AddUGrdParticle(part);
      }
      cout << "out of which " << muons_ugrd << " are muons\n";
      fSteppingAction->fPartColl.clear();

    }
  }

  return eSuccess;
}


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

  // Get the MEvent
  if (!event.HasMEvent())
    event.MakeMEvent();
  MEvent& mEvent = event.GetMEvent();

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

  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 G4XTankConstruction(fDepth); // reads fCurrentDetectorStation

    // Construct detector geometry and register w/ RunManager
    fRunManager->SetUserInitialization(fgTankConstruction); // reads fCurrentDetectorStation

    fRunManager->Initialize();
    fDetectorConstructed = true;
  }
  fgPMTAction = new G4XTankPMT();

  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 << "On ground -> StatId " << stationId << ':' << numParticles << " particle(s)" << flush;
    else
      continue;

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

    ConstructTraces(*sIt);

    // Check what happened in the associated muon counter
    if (!fCurrentDetectorStation->ExistsAssociatedCounter())
      continue;

    // Get the associated muon counter ID
    const int mId = fCurrentDetectorStation->GetAssociatedCounterId();

    // Get the counter
    if (!mEvent.HasCounter(mId))
      mEvent.MakeCounter(mId);

    mevt::Counter& eCounter = mEvent.GetCounter(mId);

    // Get the counterSimData
    if (!eCounter.HasSimData())
      eCounter.MakeSimData();

    CounterSimData& cSim = eCounter.GetSimData();

    const CoordinateSystemPtr dStationCS =
      sDetector.GetStation(stationId).GetLocalCoordinateSystem();

    // Get the muon counter depth
    fSteppingAction->fCounterDepth = -fDepth;

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

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

      fCurrentParticleIt = pIt;
      smartPointer parent(new Particle(*fCurrentParticleIt));
      fSteppingAction->SetImpingingParticle(parent);

      fgPMTAction->SetCurrentTank(fCurrentEventStationIt);
      fgPMTAction->SetParticle(*fCurrentParticleIt);

      if (fStoreOnGrdPart) {
        Particle part = *fCurrentParticleIt;
        cSim.AddGrdParticle(part);
      }

      fRunManager->BeamOn(1);

    } // particles

    cout << fSteppingAction->fPartColl.size()  << " particles underground at the counter's depth "
         << fSteppingAction->fCounterDepth/CLHEP::m << " m"
            " on station number " << sIt->GetId() << " ";

    const CoordinateSystemPtr CounterCS = fCurrentDetectorStation->GetLocalCoordinateSystem();
    G4XTankSteppingAction::PartCollectionIt pCollIt;

    unsigned int muons_ugrd = 0;
    // Fill the CounterSimData with the particles which reached the counters depth
    for (pCollIt = fSteppingAction->fPartColl.begin(); pCollIt != fSteppingAction->fPartColl.end(); ++pCollIt) {
      PartData pData = (*pCollIt);
      /*
      cout << "pData.Pname "   << pData.Pname   << " - " << pData.KinE/CLHEP::MeV << " (MeV)" << endl;
      cout << "pData.PDGCode " << pData.PDGCode << " - " << pData.Time/CLHEP::ns  << " (ns)" << endl;
      cout << "pData.PParent->GetType() " << pData.PParent->GetType() << " - " << (pData.PParent->GetType()==utl::Particle::eMuon||pData.PParent->GetType()==utl::Particle::eAntiMuon?"SHOWER MUON!!!!":"OTHER PARTICLE")<< endl;
      cout << "pData.TrackID " << pData.TrackID << endl;

      cout << "pData.x " << pData.x/CLHEP::m << endl;
      cout << "pData.y " << pData.y/CLHEP::m << endl;
      cout << "pData.z " << pData.z/CLHEP::m << endl;

      cout << "pData.px " << pData.px << endl;
      cout << "pData.py " << pData.py << endl;
      cout << "pData.pz " << pData.pz << endl;
      */
      if (pData.PParent->GetType() == utl::Particle::eMuon ||
          pData.PParent->GetType() == utl::Particle::eAntiMuon)
        ++muons_ugrd;

      const Point pPos(pData.x / CLHEP::m * utl::m, pData.y / CLHEP::m * utl::m, pData.z / CLHEP::m * utl::m, CounterCS);
      Vector pDir(pData.px, pData.py, pData.pz, CounterCS);
      // Creates new particle
      Particle part(pData.PDGCode, utl::Particle::eShower, pPos, pDir, pData.Time / CLHEP::ns * utl::ns,
                    1, pData.KinE/CLHEP::MeV * utl::MeV );
      part.SetParent(*pData.PParent);
      cSim.AddUGrdParticle(part);
    }

    cout << "out of which " << muons_ugrd << " are muons\n";
    fSteppingAction->fPartColl.clear();

  } // stations

  G4cerr << endl;

  delete fgPMTAction;
  fgPMTAction = nullptr;

  return eSuccess;
}


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


void
G4XTankSimulator::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();
    }
}