GGSLostEnergyAction.cpp

Go to the documentation of this file.

/*
 * GGSLostEnergyAction.cpp
 *
 *  Created on: 16 Jun 2011
 *      Author: Nicola mori
 */

#include "montecarlo/GGSRunManager.h"
#include "montecarlo/useractions/GGSLostEnergyAction.h"
#include "montecarlo/dataobjs/GGSTHadrIntInfo.h"
#include "montecarlo/useractions/GGSLostEnergyMessenger.h"
#include "montecarlo/GGSDetectorConstruction.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"

#include "G4StepPoint.hh"
#include "G4Box.hh"
#include "G4SystemOfUnits.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

RegisterUA(GGSLostEnergyAction);

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

GGSLostEnergyAction::GGSLostEnergyAction() :
    _outBase(""), _outTreeName(""), _outRootFile(NULL), _outTree(NULL), _kinetic(true), _tolerance(0.1 * cm) {

  _messenger = new GGSLostEnergyMessenger(this);
  _lostEnergyInfo = new GGSTLostEnergyInfo;

  GGSDetectorConstruction *detConstruction =
      (GGSDetectorConstruction*) (GGSRunManager::GetRunManager()->GetUserDetectorConstruction());
  G4ThreeVector physicalWorldPos = detConstruction->GetPhysicalWorld()->GetObjectTranslation();
  G4Box *solidWorld = (G4Box*) (detConstruction->GetPhysicalWorld()->GetLogicalVolume()->GetSolid());

  _worldXMin = physicalWorldPos[0] - solidWorld->GetXHalfLength();
  _worldXMax = physicalWorldPos[0] + solidWorld->GetXHalfLength();
  _worldYMin = physicalWorldPos[1] - solidWorld->GetYHalfLength();
  _worldYMax = physicalWorldPos[1] + solidWorld->GetYHalfLength();
  _worldZMin = physicalWorldPos[2] - solidWorld->GetZHalfLength();
  _worldZMax = physicalWorldPos[2] + solidWorld->GetZHalfLength();

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

GGSLostEnergyAction::~GGSLostEnergyAction() {

  delete _messenger;
  delete _lostEnergyInfo;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void GGSLostEnergyAction::PostUserTrackingAction(const G4Track *track) {

  const G4StepPoint *stepPoint = track->GetStep()->GetPostStepPoint();

  // Check only tracks ending because of transportation
  if (stepPoint->GetProcessDefinedStep()->GetProcessType() == fTransportation) {
    const G4ThreeVector position = stepPoint->GetPosition();
    bool charged = track->GetParticleDefinition()->GetPDGCharge() != 0;
    Double_t energy;
    if (_kinetic) {
      energy = stepPoint->GetKineticEnergy();
      _lostEnergyInfo->energyType = 1;
    }
    else {
      energy = stepPoint->GetTotalEnergy();
      _lostEnergyInfo->energyType = 0;
    }

    // Check the escape side. The sequence of checks is done to minimize the checks
    // (ie., supposing that particles are shot towards negative Z the most likely leakage
    //  is from bottom, so bottom is checked first and top last).
    // Bottom
    if (fabs(position[2] - _worldZMin) < _tolerance)
      if (charged) {
        _lostEnergyInfo->chargedBottom += energy;
        _lostEnergyInfo->nCharged++;
      }
      else {
        _lostEnergyInfo->neutralBottom += energy;
        _lostEnergyInfo->nNeutrals++;
      }
    // Front
    else if (fabs(position[0] - _worldXMax) < _tolerance)
      if (charged) {
        _lostEnergyInfo->chargedFront += energy;
        _lostEnergyInfo->nCharged++;
      }
      else {
        _lostEnergyInfo->neutralFront += energy;
        _lostEnergyInfo->nNeutrals++;
      }
    // Back
    else if (fabs(position[0] - _worldXMin) < _tolerance)
      if (charged) {
        _lostEnergyInfo->chargedBack += energy;
        _lostEnergyInfo->nCharged++;
      }
      else {
        _lostEnergyInfo->neutralBack += energy;
        _lostEnergyInfo->nNeutrals++;
      }
    // Right
    else if (fabs(position[1] - _worldYMax) < _tolerance)
      if (charged) {
        _lostEnergyInfo->chargedRight += energy;
        _lostEnergyInfo->nCharged++;
      }
      else {
        _lostEnergyInfo->neutralRight += energy;
        _lostEnergyInfo->nNeutrals++;
      }
    // Left
    else if (fabs(position[1] - _worldYMin) < _tolerance)
      if (charged) {
        _lostEnergyInfo->chargedLeft += energy;
        _lostEnergyInfo->nCharged++;
      }
      else {
        _lostEnergyInfo->neutralLeft += energy;
        _lostEnergyInfo->nNeutrals++;
      }
    // Top
    else if (fabs(position[2] - _worldZMax) < _tolerance) {
      if (charged) {
        _lostEnergyInfo->chargedTop += energy;
        _lostEnergyInfo->nCharged++;
      }
      else {
        _lostEnergyInfo->neutralTop += energy;
        _lostEnergyInfo->nNeutrals++;
      }
    }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void GGSLostEnergyAction::BeginOfEventAction(const G4Event *) {

  _lostEnergyInfo->Reset();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void GGSLostEnergyAction::EndOfEventAction(const G4Event *) {

  // Convert energies to GeV
  _lostEnergyInfo->chargedTop /= GeV;
  _lostEnergyInfo->chargedBottom /= GeV;
  _lostEnergyInfo->chargedLeft /= GeV;
  _lostEnergyInfo->chargedRight /= GeV;
  _lostEnergyInfo->chargedFront /= GeV;
  _lostEnergyInfo->chargedBack /= GeV;
  _lostEnergyInfo->neutralTop /= GeV;
  _lostEnergyInfo->neutralBottom /= GeV;
  _lostEnergyInfo->neutralLeft /= GeV;
  _lostEnergyInfo->neutralRight /= GeV;
  _lostEnergyInfo->neutralFront /= GeV;
  _lostEnergyInfo->neutralBack /= GeV;

  if (_outTreeName != "") {
    // Call fill only if not using the default tree
    _outRootFile->cd();
    _outTree->Fill();
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void GGSLostEnergyAction::BeginOfRunAction(const G4Run *run) {

  _outRootFile = GGSRootFileService::GetInstance().GetFileForThisRun(_outBase, run);
  if (_outTreeName == "") {
    _outTree = GGSRootFileService::GetInstance().GetDefaultTree(_outRootFile);
    _outTree->SetTitle(TString(_outTree->GetTitle()) + "LostEn ");
  }
  else
    _outTree = new TTree(_outTreeName.data(), "Energy lost by escaping particles");

  _outTree->Branch("lostEnergyInfo", "GGSTLostEnergyInfo", &_lostEnergyInfo);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void GGSLostEnergyAction::EndOfRunAction(const G4Run *) {

  if (_outTreeName != "") {
    // Write the tree if we are not using the default one
    _outRootFile->cd();
    _outTree->Write();
  }
  GGSRootFileService::GetInstance().CloseFileForThisRun(_outBase);
  _outRootFile = NULL;
  _outTree = NULL;
}