FDsimG4EventAction.cc

Go to the documentation of this file.

#include <fwk/RunController.h>
#include <fwk/VModule.h>
#include <fwk/RandomEngineRegistry.h>
#include <utl/RandomEngine.h>
#include "FDsimG4Defs.hh"
#include "FDsimG4EventAction.hh"
#include "FDsimG4DetectorConstruction.hh"
#include "FDsimG4PrimaryGeneratorAction.hh"
#include "FDsimG4Write2ROOT.hh"
#include "G4RunManager.hh"
#include "G4Run.hh"
#include "G4Event.hh"
#include "G4EventManager.hh"
#include "G4SDManager.hh"
#include "G4GeneralParticleSource.hh"
#include "G4TrajectoryContainer.hh"
#include "G4Trajectory.hh"
#include "G4VVisManager.hh"
#include <cmath>
#include <fstream>

using namespace std;
using namespace TelescopeSimulatorLX;
using CLHEP::RandFlat;


FDsimG4EventAction::FDsimG4EventAction()
{
  Init();
}

  
FDsimG4EventAction::FDsimG4EventAction(FDsimG4Write2ROOT* write2Root) 
{
  Init();
  fWrite2ROOT = write2Root;
}


FDsimG4EventAction::~FDsimG4EventAction()
{ }


void
FDsimG4EventAction::Init()
{
  fOpticalHitsCollID = -1;
  fOpticalHitsColl = NULL;
  fFDSimulator = NULL;
  fWrite2ROOT = NULL; 

  TelescopeSimulator& FDSimulator =
    dynamic_cast<TelescopeSimulator&>(fwk::RunController::GetInstance().GetModule("TelescopeSimulatorLX"));
  fFDSimulator = &FDSimulator;
}


void
FDsimG4EventAction::BeginOfEventAction(const G4Event* const)
{
#ifdef __TELESCOPE_SIMULATOR_LX_USE_ROOT__
  if (fWrite2ROOT) 
    fWrite2ROOT->Clean();
#endif
}


void
FDsimG4EventAction::EndOfEventAction(const G4Event* const evt)
{
  static G4SDManager* SDman = G4SDManager::GetSDMpointer();

  if (fOpticalHitsCollID == -1)
    fOpticalHitsCollID = SDman->GetCollectionID("OpticalHits");

  G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
  if (HCE) {
    if (fOpticalHitsCollID != -1) 
      fOpticalHitsColl = (FDsimG4OpticalHitsCollection*)(HCE->GetHC(fOpticalHitsCollID));
  }

  G4String drawFlag = "all";
  G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
  if (pVVisManager) {
    G4TrajectoryContainer* trajectoryContainer = evt->GetTrajectoryContainer();
    G4int n_trajectories = 0;
    if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();
    for (G4int i = 0; i < n_trajectories; i++){ 
      G4Trajectory* trj = (G4Trajectory *)((*(evt->GetTrajectoryContainer()))[i]);
      if (drawFlag == "all") 
        trj->DrawTrajectory();
      else if ((drawFlag == "charged") && (trj->GetCharge() != 0.))
        trj->DrawTrajectory();
    }
  }

  if (fOpticalHitsColl) 
    fOpticalHitsColl->DrawAllHits();

  if (fFDSimulator && fOpticalHitsColl) {
      for (G4int iHit = 0; iHit < fOpticalHitsColl->entries(); iHit++){
        FDsimG4StoreOpticalHit StoreOpticalHit = FDsimG4StoreOpticalHit();
        CopyHit((*fOpticalHitsColl)[iHit], &StoreOpticalHit);
        fFDSimulator->SetG4PixelHit(StoreOpticalHit);
      }
  }

#ifdef __TELESCOPE_SIMULATOR_LX_USE_ROOT__
  static G4bool useSD =
    (dynamic_cast<const FDsimG4DetectorConstruction*>
     ((G4RunManager::GetRunManager())->GetUserDetectorConstruction()))->
    GetUseSensitiveDetectors();

  if (fWrite2ROOT && useSD && HCE){
    WriteHits("FilterHits",HCE);
    WriteHits("LensHits",HCE);
    WriteHits("MirrorHits",HCE);
    WriteHits("MercedesHits",HCE);
    WriteHits("FocalSurfaceHits",HCE);
    WriteHits("OpticalHits",HCE);

    FDsimG4StoreOpticalHit thePrimary = FDsimG4StoreOpticalHit();

    G4PrimaryVertex* vertex = evt->GetPrimaryVertex();
    G4PrimaryParticle* particle = vertex->GetPrimary();

    G4double T0 = vertex->GetT0();
    G4double X0 = vertex->GetX0();
    G4double Y0 = vertex->GetY0();
    G4double Z0 = vertex->GetZ0();
    G4double weight = vertex->GetWeight();

    G4double p0 = particle->GetMomentum().mag();
    G4double mass = particle->GetMass();
    G4double E0 = sqrt(p0*p0+mass*mass);

    G4double DirX = particle->GetMomentum().x()/p0;
    G4double DirY = particle->GetMomentum().y()/p0;
    G4double DirZ = particle->GetMomentum().z()/p0;

    thePrimary.Reset();
    thePrimary.SetEnergy(E0/eV);
    thePrimary.SetTime(T0/ns);
    thePrimary.SetWavelength(h_Planck*c_light/E0/nm);
    thePrimary.SetPosX(X0/mm);
    thePrimary.SetPosY(Y0/mm);
    thePrimary.SetPosZ(Z0/mm);
    
    thePrimary.SetDirX(DirX);
    thePrimary.SetDirY(DirY);
    thePrimary.SetDirZ(DirZ);
    thePrimary.SetWeight(weight);

    fWrite2ROOT->SetINPhoton(&thePrimary);

    fWrite2ROOT->WriteIt();
  }
#endif
}


#ifdef __TELESCOPE_SIMULATOR_LX_USE_ROOT__

void
FDsimG4EventAction::WriteHits(string HitsCollectionName,G4HCofThisEvent* HCE)
{
  static G4SDManager* SDman = G4SDManager::GetSDMpointer();
  G4int HitsCollID  = SDman->GetCollectionID(HitsCollectionName);
  FDsimG4OpticalHitsCollection* HitsCollection;

  HitsCollection = (FDsimG4OpticalHitsCollection*)(HCE->GetHC(HitsCollID));
  if (!HitsCollection) 
    return;

  G4int nHits = HitsCollection->entries();
  if (nHits>0){
    G4int iHit = nHits-1;

    FDsimG4StoreOpticalHit StoreOpticalHit = FDsimG4StoreOpticalHit();

    CopyHit((*HitsCollection)[iHit], &StoreOpticalHit);
    

    if (HitsCollectionName == "FilterHits")
      fWrite2ROOT->SetG4FilterPhoton(&StoreOpticalHit);
    else
      if (HitsCollectionName == "LensHits")
        fWrite2ROOT->SetG4LensPhoton(&StoreOpticalHit);
      else
        if (HitsCollectionName == "MirrorHits")
          fWrite2ROOT->SetG4MirrorPhoton(&StoreOpticalHit);
        else
          if (HitsCollectionName == "MercedesHits")
            fWrite2ROOT->SetG4MercedesPhoton(&StoreOpticalHit);
          else
            if (HitsCollectionName == "FocalSurfaceHits")
              fWrite2ROOT->SetG4FSPhoton(&StoreOpticalHit);  
            else
              if (HitsCollectionName == "OpticalHits")
                fWrite2ROOT->SetG4Photon(&StoreOpticalHit);
  }
}
#endif


void 
FDsimG4EventAction::CopyHit(const FDsimG4OpticalHit* theHit, 
                            FDsimG4StoreOpticalHit* theStoreHit)
{
  G4double HitEnergy = theHit->GetEnergy();
  G4double HitTime = theHit->GetTime();
  G4double HitWavelength = theHit->GetWavelength();
  G4ThreeVector HitPos = theHit->GetWorldPosition();
  G4ThreeVector HitLocalPos = theHit->GetPosition();
  G4ThreeVector PixelPos = theHit->GetPixelPosition();
  G4ThreeVector HitDir = theHit->GetDirection();
  G4int HitID = theHit->GetHitID();
  G4int PMTID = theHit->GetPMTid();
  G4double weight = theHit->GetWeight();

  theStoreHit->Reset();
  theStoreHit->SetEnergy(HitEnergy/eV);
  theStoreHit->SetTime(HitTime/ns);
  theStoreHit->SetWavelength(HitWavelength/nanometer);
  theStoreHit->SetPosX(HitPos.x()/mm);
  theStoreHit->SetPosY(HitPos.y()/mm);
  theStoreHit->SetPosZ(HitPos.z()/mm);

  static utl::RandomEngine& RandomEngine = 
    fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);

  double r = HitPos.mag() + 
    0.000001*mm * RandFlat::shoot(&RandomEngine.GetEngine(),0.0,1.0);

  theStoreHit->SetPosR(r/mm);
  theStoreHit->SetDirX(HitDir.x());
  theStoreHit->SetDirY(HitDir.y());
  theStoreHit->SetDirZ(HitDir.z());
  theStoreHit->SetLocalPosX(HitLocalPos.x()/mm);
  theStoreHit->SetLocalPosY(HitLocalPos.y()/mm);
  theStoreHit->SetLocalPosZ(HitLocalPos.z()/mm);
  theStoreHit->SetPixelPosX(PixelPos.x()/mm);
  theStoreHit->SetPixelPosY(PixelPos.y()/mm);
  theStoreHit->SetPixelPosZ(PixelPos.z()/mm);
  theStoreHit->SetHitID(HitID);
  theStoreHit->SetPMTid(PMTID);
  theStoreHit->SetWeight(weight);
}