GGSMultiParticleGeneratorAction.cpp

Go to the documentation of this file.

/*
 * GGSMultiParticleGeneratorAction.cpp
 *
 *  Created on: 09/may/2011
 *      Author: Nicola mori
 */

#include "montecarlo/generators/GGSMultiParticleGeneratorAction.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"

#include "G4Event.hh"
#include "G4SystemOfUnits.hh"

RegisterGA(GGSMultiParticleGeneratorAction, multiParticle);

GGSMultiParticleGeneratorAction::GGSMultiParticleGeneratorAction()
    : _messenger(this, "/GGS/generatorActions/multiParticle/") {

  _messenger
      .DeclareMethod("file", &GGSMultiParticleGeneratorAction::SetEventsFile,
                     "Set the name of the file containing the definitions of the events")
      .SetStates(G4State_Idle);
}

GGSMultiParticleGeneratorAction::~GGSMultiParticleGeneratorAction() {
  if (_inputFile.is_open()) {
    _inputFile.close();
  }
}

void GGSMultiParticleGeneratorAction::SetEventsFile(const G4String &evFile) {
  if (_inputFile.is_open()) {
    _inputFile.close();
  }
  _inputFile.open(evFile.data());
  if (!_inputFile) {
    G4Exception("GGSMultiParticleGeneratorAction::GGSMultiParticleGeneratorAction", "mp-nofile", FatalException,
                "Cannot open file");
  }
}

// This member function is largely copied from G4HEPEvtInterface::GeneratePrimaryVertex
// Instead of relying on a G4VParticleGenerator, generate the primaries directly here
void GGSMultiParticleGeneratorAction::GeneratePrimaries(G4Event *anEvent) {

  G4int NHEP; // number of entries
  _inputFile >> NHEP;
  if (_inputFile.eof()) {
    G4Exception("GGSMultiParticleGeneratorAction::GeneratePrimaryVertex", "mp-eof", JustWarning,
                "End-Of-File: multi-particle input file");
    return;
  }

  G4int ISTHEP;   // status code
  G4int IDHEP;    // PDG code
  G4int JDAHEP1;  // first daughter
  G4int JDAHEP2;  // last daughter
  G4double PHEP1; // px in GeV
  G4double PHEP2; // py in GeV
  G4double PHEP3; // pz in GeV
  G4double PHEP5; // mass in GeV
  G4double VHEP1; // x position in cm
  G4double VHEP2; // y position in cm
  G4double VHEP3; // z position in cm
  // WARNING: in standard HEPEVT, VHEP4 is in mm/c. Here we use ns since it's more comfortable within Geant4
  G4double VHEP4; // production time in ns

  for (G4int IHEP = 0; IHEP < NHEP; IHEP++) {

    _inputFile >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2 >> PHEP1 >> PHEP2 >> PHEP3 >> PHEP5 >> VHEP1 >> VHEP2 >>
        VHEP3 >> VHEP4;

    // create G4PrimaryParticle object
    G4PrimaryParticle *particle = new G4PrimaryParticle(IDHEP, PHEP1 * GeV, PHEP2 * GeV, PHEP3 * GeV);
    particle->SetMass(PHEP5 * GeV);

    // create G4PrimaryVertex object
    G4PrimaryVertex *vertex = new G4PrimaryVertex(G4ThreeVector(VHEP1 * cm, VHEP2 * cm, VHEP3 * cm), VHEP4 * ns);

    // Create Particle object
    _particles.emplace_back(Particle());
    Particle &hepParticle = _particles.back();
    hepParticle.theParticle = particle;
    hepParticle.ISTHEP = ISTHEP;
    hepParticle.JDAHEP1 = JDAHEP1;
    hepParticle.JDAHEP2 = JDAHEP2;
    hepParticle.theVertex = vertex;
  }

  // check if there is at least one particle
  if (_particles.size() == 0)
    return;

  // make connection between daughter particles decayed from
  // the same mother
  for (size_t iPart = 0; iPart < _particles.size(); ++iPart) {
    if (_particles[iPart].JDAHEP1 > 0) //  it has daughters
    {
      G4int jda1 = _particles[iPart].JDAHEP1 - 1; // FORTRAN index starts from 1
      G4int jda2 = _particles[iPart].JDAHEP2 - 1; // but C++ starts from 0.
      G4PrimaryParticle *mother = _particles[iPart].theParticle;
      for (G4int jPart = jda1; jPart <= jda2; ++jPart) {
        G4PrimaryParticle *daughter = _particles[jPart].theParticle;
        if (_particles[jPart].ISTHEP > 0) {
          mother->SetDaughter(daughter);
          _particles[jPart].ISTHEP *= -1;
        }
      }
    }
  }

  // put initial particles to the vertex and the vertex to G4Event object
  for (size_t iPart = 0; iPart < _particles.size(); ++iPart) {
    if (_particles[iPart].ISTHEP > 0) { // ISTHEP of daughters had been set to negative
      _particles[iPart].theVertex->SetPrimary(_particles[iPart].theParticle);
      anEvent->AddPrimaryVertex(_particles[iPart].theVertex);
    }
  }

  // clear G4HEPEvtParticles
  _particles.clear();
}

GGSParameters GGSMultiParticleGeneratorAction::GetParameters() const {
  GGSParameters params;
  params.SetParam("generator", std::string("multiParticle"));
  return params;
}