GGSHadrIntAction.cpp

Go to the documentation of this file.

/*
 * GGSHadrIntAction.cpp
 *
 *  Created on: 31 May 2011
 *      Author: Nicola mori
 */

#include "montecarlo/useractions/GGSHadrIntAction.h"
#include "montecarlo/dataobjs/GGSTHadrIntInfo.h"
#include "montecarlo/dataobjs/GGSTParticle.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"
#include "montecarlo/useractions/GGSHadrIntMessenger.h"

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

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

RegisterUA(GGSHadrIntAction);

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

GGSHadrIntAction::GGSHadrIntAction()
    : _energyFraction(0.9), _outputProducts(false), _outBase(""), _outTreeName(""), _outRootFile(NULL), _outTree(NULL) {

  _messenger = new GGSHadrIntMessenger(this);
  _quasiElasticInteractions = new TClonesArray("GGSTHadrIntInfo", 10);
  _inelasticInteractions = new TClonesArray("GGSTHadrIntInfo", 1);
}

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

GGSHadrIntAction::~GGSHadrIntAction() {

  delete _messenger;
  delete _quasiElasticInteractions;
  delete _inelasticInteractions;
}

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

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

  // Check whether it is a primary particle
  G4PrimaryParticle *primary = track->GetDynamicParticle()->GetPrimaryParticle();
  if (primary != NULL)
    // Check whether it is a baryon
    if (primary->GetParticleDefinition()->GetBaryonNumber() > 0)
      _tracks[track] = track->GetTrackID();
}

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

void GGSHadrIntAction::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()) {
      // Check process type
      const G4Step *step = track->GetStep();
      int pType = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessType();
      if (pType == fHadronic) {
        int primaryPDG = track->GetDefinition()->GetPDGEncoding();

        // Check if interaction is quasi-elastic
        float maxSecondaryEnergy = step->GetPreStepPoint()->GetKineticEnergy() * _energyFraction;
        auto *secondaries = step->GetSecondaryInCurrentStep();
        unsigned int nSecondaries = secondaries->size();
        bool inelastic = true;
        GGSTHadrIntInfo *newElement = NULL;
        for (unsigned int is = 0; is < nSecondaries; is++) {
          auto *currentSecondaryTrack = (*secondaries)[is];
          if (currentSecondaryTrack->GetDefinition()->GetPDGEncoding() == primaryPDG &&
              currentSecondaryTrack->GetKineticEnergy() > maxSecondaryEnergy) {

            // Quasi-elastic interaction
            newElement =
                (GGSTHadrIntInfo *)(_quasiElasticInteractions->ConstructedAt(_quasiElasticInteractions->GetEntries()));
            newElement->Clear();

            // Add an entry in the track list for the new "primary"
            _tracks[currentSecondaryTrack] = trackIter->second;
            // Set flag to avoid treating it as inelastic
            inelastic = false;
            // End loop
            break;
          }
        }

        if (inelastic) {
          // Inelastic interaction
          newElement = (GGSTHadrIntInfo *)(_inelasticInteractions->ConstructedAt(_inelasticInteractions->GetEntries()));
          newElement->Clear();
        }

        // Fill the new element
        newElement->primary.PDGCode = primaryPDG;
        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;
        const G4Isotope *target =
            ((G4HadronicProcess *)(step->GetPostStepPoint()->GetProcessDefinedStep()))->GetTargetIsotope();
        newElement->targetPDGCode = 1000000000 + target->GetZ() * 10000 + target->GetN() * 10;
        newElement->originalTrackID = trackIter->second;
        newElement->processName = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName().data();
        newElement->volumeName = step->GetPostStepPoint()->GetPhysicalVolume()->GetName();
        if (_outputProducts) {
          if (newElement->products) {
            newElement->products->Clear("C");
          } else {
            newElement->products = new TClonesArray("GGSTParticle");
          }
          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 = NULL;
        } // if(_outputProducts)
      }   // if (pType == fHadronic)

      // Delete the currently followed particle from list
      _tracks.erase(trackIter);
    } // if (trackIter != _tracks.end())
  }
}

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

void GGSHadrIntAction::BeginOfEventAction(const G4Event *) {

  _inelasticInteractions->Clear("C");
  _quasiElasticInteractions->Clear("C");
}

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

void GGSHadrIntAction::EndOfEventAction(const G4Event *) {

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

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

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

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

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

  // Create quasi-elastic interactions branch
  _outTree->Branch("quasiElastic", "TClonesArray", &_quasiElasticInteractions);
}

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

void GGSHadrIntAction::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;
}