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 = new G4GenericMessenger(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();
  }
  delete _messenger;

}

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
    Particle* hepParticle = new Particle;
    hepParticle->theParticle = particle;
    hepParticle->ISTHEP = ISTHEP;
    hepParticle->JDAHEP1 = JDAHEP1;
    hepParticle->JDAHEP2 = JDAHEP2;
    hepParticle->theVertex = vertex;

    // Store
    _particles.push_back(hepParticle);
  }

// 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 i = 0; i < _particles.size(); i++) {
    if (_particles[i]->JDAHEP1 > 0) //  it has daughters
        {
      G4int jda1 = _particles[i]->JDAHEP1 - 1; // FORTRAN index starts from 1
      G4int jda2 = _particles[i]->JDAHEP2 - 1; // but C++ starts from 0.
      G4PrimaryParticle* mother = _particles[i]->theParticle;
      for (G4int j = jda1; j <= jda2; j++) {
        G4PrimaryParticle* daughter = _particles[j]->theParticle;
        if (_particles[j]->ISTHEP > 0) {
          mother->SetDaughter(daughter);
          _particles[j]->ISTHEP *= -1;
          ;
        }
      }
    }
  }

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

// clear G4HEPEvtParticles
//_particles.clearAndDestroy();
  for (size_t iii = 0; iii < _particles.size(); iii++) {
    delete _particles[iii];
  }
  _particles.clear();

//evt->AddPrimaryVertex(vertex);

}