GGSPrimaryDisAction.cpp

Go to the documentation of this file.

/*
 * GGSPrimaryDisAction.cpp
 *
 *  Created on: 1 Jun 2020
 *      Author: Nicola Mori
 */

#include "montecarlo/useractions/GGSPrimaryDisAction.h"
#include "montecarlo/dataobjs/GGSTParticle.h"
#include "montecarlo/dataobjs/GGSTPrimaryDisInfo.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"

#include "G4HadronicProcess.hh"
#include "G4SystemOfUnits.hh"

RegisterUA(GGSPrimaryDisAction);

GGSPrimaryDisAction::GGSPrimaryDisAction()
    : _disInfo("GGSTPrimaryDisInfo", 1), _outputProducts(false), _outBase(""), _outTreeName(""), _outRootFile(nullptr),
      _outTree(nullptr),
      _messenger(this, "/GGS/userActions/primaryDisAction/", "Commands for primary disappearance action") {

  _messenger.DeclareProperty("fileBase", _outBase, "Sets the base name for ROOT output file")
      .SetGuidance("   Can be with or without extension (.root will be used automatically)")
      .SetStates(G4State_PreInit, G4State_Idle);
  _messenger.DeclareProperty("treeName", _outTreeName, "Set the name of the TTree object in the ROOT output file.")
      .SetStates(G4State_PreInit, G4State_Idle);
  _messenger.DeclareProperty("outProducts", _outputProducts, "Toggles on and off saving interaction products.")
      .SetParameterName("outProducts", false)
      .SetStates(G4State_PreInit, G4State_Idle);
}

void GGSPrimaryDisAction::PreUserTrackingAction(const G4Track *track) {

  // Check whether it is a primary particle
  G4PrimaryParticle *primary = track->GetDynamicParticle()->GetPrimaryParticle();
  if (primary != nullptr) {
    _tracks[track] = track->GetTrackID();
  }
}

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

  if (_tracks.size() > 0) {
    // See if the track is in the tracking list
    std::map<const G4Track *, int>::iterator trackIter = _tracks.find(track);
    if (trackIter != _tracks.end()) {
      GGSTPrimaryDisInfo *newElement = (GGSTPrimaryDisInfo *)(_disInfo.ConstructedAt(_disInfo.GetEntries()));
      newElement->Clear();

      // Fill the new element
      newElement->primary.PDGCode = track->GetDefinition()->GetPDGEncoding();
      const G4Step *step = track->GetStep();
      newElement->primary.trackID = step->GetTrack()->GetTrackID();
      for (int i = 0; i < 3; i++) {
        newElement->primary.pos[i] = track->GetPosition()[i] / cm;
        newElement->primary.mom[i] = track->GetMomentum()[i] / GeV;
      }
      newElement->primary.time = track->GetGlobalTime() / ns;
      newElement->processName = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName().data();
      G4VPhysicalVolume *postStepVol = step->GetPostStepPoint()->GetPhysicalVolume();
      if (postStepVol) {
        newElement->volumeName = postStepVol->GetName();
      } else {
        newElement->volumeName = "outOfWorld";
      }
      if (_outputProducts) {
        if (newElement->products) {
          newElement->products->Clear("C");
        } else {
          newElement->products = new TClonesArray("GGSTParticle");
        }
        auto *secondaries = step->GetSecondaryInCurrentStep();
        unsigned int nSecondaries = secondaries->size();
        for (unsigned int is = 0; is < nSecondaries; is++) {
          GGSTParticle *sec = new ((*(newElement->products))[is]) GGSTParticle;
          auto *currentSecondaryTrack = (*secondaries)[is];
          sec->PDGCode = currentSecondaryTrack->GetDefinition()->GetPDGEncoding();
          sec->trackID = currentSecondaryTrack->GetTrackID();
          for (int i = 0; i < 3; i++) {
            sec->pos[i] = currentSecondaryTrack->GetPosition()[i] / cm;
            sec->mom[i] = currentSecondaryTrack->GetMomentum()[i] / GeV;
          }
          sec->time = currentSecondaryTrack->GetGlobalTime() / ns;
        }

      } else {
        delete newElement->products;
        newElement->products = nullptr;
      } // if(_outputProducts)

      // Delete the currently followed particle from list
      _tracks.erase(trackIter);

    } // if (trackIter != _tracks.end())

  } // if (_tracks.size() > 0)
}

void GGSPrimaryDisAction::BeginOfEventAction(const G4Event *) { _disInfo.Clear("C"); }

void GGSPrimaryDisAction::EndOfEventAction(const G4Event *) {

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

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

  _outRootFile = GGSRootFileService::GetInstance().GetFileForThisRun(_outBase, run);
  if (_outTreeName == "") {
    _outTree = GGSRootFileService::GetInstance().GetDefaultTree(_outRootFile);
    _outTree->SetTitle(TString(_outTree->GetTitle()) + "DisInfo ");
  } else
    _outTree = new TTree(_outTreeName.data(), "Disappearance of primary particles");

  // Create inelastic interactions branch
  _outTree->Branch("primaryDisInfo", "TClonesArray", &_disInfo);
}

void GGSPrimaryDisAction::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 = nullptr;
  _outTree = nullptr;
}