GGSMCTruthAction.cpp

Go to the documentation of this file.

/*
 * GGSMCTruthAction.cpp
 *
 *  Created on: 17 Aug 2011
 *      Author: Nicola mori
 */

// GGS headers
#include "montecarlo/useractions/GGSMCTruthAction.h"
#include "montecarlo/GGSRunManager.h"
#include "montecarlo/dataobjs/GGSPrimaryParticleInfo.h"
#include "montecarlo/dataobjs/GGSTMCTruthInfo.h"
#include "montecarlo/generators/GGSGunGeneratorAction.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"
#include "montecarlo/useractions/GGSMCTruthMessenger.h"

#include "G4SystemOfUnits.hh"

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

RegisterUA(GGSMCTruthAction);

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

GGSMCTruthAction::GGSMCTruthAction() : _outBase(""), _outTreeName(""), _outRootFile(NULL), _outTree(NULL) {

  _messenger = new GGSMCTruthMessenger(this);
  _mcTruthInfo = new GGSTMCTruthInfo;
  _primaryParticle = new GGSTParticle;
}

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

GGSMCTruthAction::~GGSMCTruthAction() {

  delete _messenger;
  delete _mcTruthInfo;
  delete _primaryParticle;
}

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

void GGSMCTruthAction::EndOfEventAction(const G4Event *event) {

  _Convert(event, _mcTruthInfo);

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

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

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

  _outRootFile = GGSRootFileService::GetInstance().GetFileForThisRun(_outBase, run);
  if (_outTreeName == "") {
    _outTree = GGSRootFileService::GetInstance().GetDefaultTree(_outRootFile);
    _outTree->SetTitle(TString(_outTree->GetTitle()) + "MCTruth ");
  } else
    _outTree = new TTree(_outTreeName.data(), "MC truth");

  // Create MC truth info branch
  _outTree->Branch("mcTruthInfo", "GGSTMCTruthInfo", &_mcTruthInfo);
}

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

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

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

void GGSMCTruthAction::_Convert(const G4Event *event, GGSTMCTruthInfo *&mcTruth) {

  // -----------------------------------------------
  // check if GGSTMCTruthInfo exists. if not, create it
  // -----------------------------------------------
  if (!mcTruth)
    mcTruth = new GGSTMCTruthInfo();
  else
    mcTruth->Reset();

  mcTruth->eventID = event->GetEventID();

  // Obtain the number of discarded events
  mcTruth->nDiscarded = GGSRunManager::GetRunManager()->GetNDiscardedEvents();
  // Add killed events to discarded ones
  mcTruth->nDiscarded += GGSRunManager::GetRunManager()->GetNKilledEvents();

  G4int n_vertex = event->GetNumberOfPrimaryVertex();
  //-----------------------
  // Loop over all vertexes
  //-----------------------
  for (G4int iv = 0; iv < n_vertex; iv++) {
    G4PrimaryVertex *pv = event->GetPrimaryVertex(iv);
    //-----------------------------------
    // Loop over primaries in this vertex
    //-----------------------------------
    for (G4int ipp = 0; ipp < pv->GetNumberOfParticle(); ipp++) {

      G4PrimaryParticle *pp = pv->GetPrimary(ipp);

      GGSTParticle *tPP = new GGSTParticle();
      tPP->pos[0] = pv->GetX0() / cm;
      tPP->pos[1] = pv->GetY0() / cm;
      tPP->pos[2] = pv->GetZ0() / cm;
      tPP->time = pv->GetT0() / ns;

      //------------------------------
      // Extract primary particle info
      //------------------------------

      tPP->PDGCode = pp->GetPDGcode();
      tPP->trackID = pp->GetTrackID();
      tPP->mom[0] = pp->GetPx() / GeV;
      tPP->mom[1] = pp->GetPy() / GeV;
      tPP->mom[2] = pp->GetPz() / GeV;

      GGSPrimaryParticleInfo *partInfo = (GGSPrimaryParticleInfo *)(pp->GetUserInformation());
      if (partInfo) {
        tPP->isTrackKilled = partInfo->isTrackKilled;
      } else {
        tPP->isTrackKilled = 2;
      }

      //-----------------------------
      // Record primary particle info
      //-----------------------------
      mcTruth->AddPrimaryParticle(tPP);
      delete tPP;
    }
  }
}

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