GGSGunGeneratorAction.cpp

Go to the documentation of this file.

/*
 * GGSGunGeneratorAction.cpp
 *
 *  Created on: 09 Oct 2013
 *      Author: Nicola Mori
 */

// GGS headers
#include "montecarlo/generators/GGSGunGeneratorAction.h"
#include "montecarlo/GGSDetectorConstruction.h"
#include "montecarlo/GGSRunManager.h"
#include "montecarlo/generators/GGSGunGeneratorActionMessenger.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"
#include "utils/GGSSmartLog.h"

// Geant4/CLHEP headers
#include "G4Event.hh"
#include "G4IonTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleGun.hh"
#include "G4SystemOfUnits.hh"
#include "G4UImanager.hh"
#include "Randomize.hh"

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

RegisterGA(GGSGunGeneratorAction, gun);

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

GGSGunGeneratorAction::GGSGunGeneratorAction()
    : _shootAxis(0, 0, 1), _acceptanceFlag(false), _rndmMinPositionFlag(false), _rndmMaxPositionFlag(false),
      _rndmFromFlat(false), _rndmFromSphere(false), _rndmSphereCenter(0., 0., 0.), _rndmSphereCapExtRand(0.),
      _rndmSphereCapAlpha(0.), _rndmSphereCapBeta(0.), _rndmSphereCapRotMat(G4RotationMatrix::IDENTITY),
      _rndmSphereRadius(0.), _rndmMinEnergyFlag(false), _rndmMaxEnergyFlag(false), _spectralIndexFlag(false),
      _rndmMinThetaFlag(false), _rndmMaxThetaFlag(false), _rndmMinPhiFlag(false), _rndmMaxPhiFlag(false),
      _rndmPosFlag(false), _rndmDirFlag(false), _position(0., 0., 0.), _minPosition(0., 0., 0.),
      _maxPosition(0., 0., 0.), _energy(10. * GeV), _minEnergy(10. * GeV), _maxEnergy(10. * GeV), _spectralIndex(-3.2),
      _sinTheta(0.), _cosTheta(1.), _minThetaRand(0.), _maxThetaRand(.5), _minCos2ThetaRand(0.), _maxCos2ThetaRand(1.),
      _sinPhi(0.), _cosPhi(1.), _minPhi(0.), _maxPhi(2 * CLHEP::pi), _direction(0., 0., 0.) {

  _SetShootingDir();

  _gunGenerator = new G4ParticleGun(1);
  _gunGenerator->SetParticleDefinition(G4ParticleTable::GetParticleTable()->FindParticle("geantino"));

  // Create the messenger for this class
  _messenger = new GGSGunGeneratorActionMessenger(this);
}

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

GGSGunGeneratorAction::~GGSGunGeneratorAction() {
  delete _gunGenerator;
  delete _messenger;
}

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

void GGSGunGeneratorAction::SetGunParticle(const G4String &particle) {
  if (particle == "C12") {
    // G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
    G4ParticleDefinition *c12 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(6, 12, 0.);
    ((G4ParticleGun *)_gunGenerator)->SetParticleDefinition(c12);
    ((G4ParticleGun *)_gunGenerator)->SetParticleCharge(6 * eplus);
  } else if (particle == "Fe56") {
    // G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
    G4ParticleDefinition *fe56 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(26, 56, 0.);
    ((G4ParticleGun *)_gunGenerator)->SetParticleDefinition(fe56);
    ((G4ParticleGun *)_gunGenerator)->SetParticleCharge(26 * eplus);
  } else
    G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gun/particle ").append(particle));
}

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

G4RotationMatrix &GGSGunGeneratorAction::_ComputeRotMatrix(double cosAlpha, double sinAlpha, double cosBeta,
                                                           double sinBeta) {

  // A small helper class which allows to directly specify the rotation matrix elements
  class Rot3 : public G4RotationMatrix {
  public:
    void SetMatrix(double mxx, double mxy, double mxz, double myx, double myy, double myz, double mzx, double mzy,
                   double mzz) {
      rxx = mxx;
      rxy = mxy;
      rxz = mxz;
      ryx = myx;
      ryy = myy;
      ryz = myz;
      rzx = mzx;
      rzy = mzy;
      rzz = mzz;
    }
  };

  static Rot3 retMatrix;

  /* The rotation is computed as the composition of the following operations (in this order):
   *
   * 1) rotate the vector so that beta = 0
   * 2) rotate about the new y axis of an angle alpha
   * 3) invert 1
   * The matrices are:
   *
   *      | cosb   sinb   0|
   *  1 = |-sinb   cosb   0|
   *      |  0      0     1|
   *
   *      | cosa    0  sina|
   *  2 = |  0      1     0|
   *      |-sina    0  cosa|
   *
   *      | cosb  -sinb   0|
   *  3 = | sinb   cosb   0|
   *      |  0      0     1|
   *
   *
   * The resulting matrix is obtained as M = 3*2*1.
   */

  double xx = cosAlpha * cosBeta * cosBeta + sinBeta * sinBeta;
  double xy = cosBeta * cosAlpha * sinBeta - cosBeta * sinBeta; // 1st row 2nd column
  double xz = cosBeta * sinAlpha;
  double yx = sinBeta * cosAlpha * cosBeta - sinBeta * cosBeta; // 2nd row 1st column
  double yy = sinBeta * cosAlpha * sinBeta + cosBeta * cosBeta;
  double yz = sinBeta * sinAlpha;
  double zx = -sinAlpha * cosBeta;
  double zy = -sinAlpha * sinBeta;
  double zz = cosAlpha;

  retMatrix.SetMatrix(xx, xy, xz, yx, yy, yz, zx, zy, zz);

  return retMatrix;
}

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

void GGSGunGeneratorAction::_ResetFlat() {

  _rndmFromFlat = false;
  _minPosition = _maxPosition = G4ThreeVector(0., 0., 0.);
  _shootAxis = G4ThreeVector(0, 0, 1);
}

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

void GGSGunGeneratorAction::_ResetSphere() {

  _rndmFromSphere = false;
  _rndmSphereCenter = G4ThreeVector(0., 0., 0.);
  _rndmSphereRadius = 0;
  _rndmSphereCapExtRand = 0.;
  _rndmSphereCapRotMat = G4RotationMatrix::IDENTITY;
}

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

void GGSGunGeneratorAction::_SetRndmFromSurface(int surfType) {

  if (surfType == 0) { // No generation from surface
    _ResetFlat();
    _ResetSphere();
  } else if (surfType == 1) { // Flat surface
    if (_minPosition[0] < _maxPosition[0] && _minPosition[1] < _maxPosition[1] && _minPosition[2] == _maxPosition[2]) {
      _rndmFromFlat = true;
      _shootAxis = G4ThreeVector(0, 0, 1);
    } else if (_minPosition[0] == _maxPosition[0] && _minPosition[1] < _maxPosition[1] &&
               _minPosition[2] < _maxPosition[2]) {
      _rndmFromFlat = true;
      _shootAxis = G4ThreeVector(1, 0, 0);
    } else if (_minPosition[0] < _maxPosition[0] && _minPosition[1] == _maxPosition[1] &&
               _minPosition[2] < _maxPosition[2]) {
      _rndmFromFlat = true;
      _shootAxis = G4ThreeVector(0, 1, 0);
    } else {
      _rndmFromFlat = false;
      _shootAxis = G4ThreeVector(0, 0, 1); // Default shoot axis along Z
    }
    if (_rndmFromFlat)
      _rndmPosFlag = true;
    _ResetSphere();
  } else if (surfType == 2) { // Spherical surface
    if (_rndmSphereRadius > 0. && _rndmSphereCapExtRand > 0.) {
      _rndmFromSphere = true;
      _rndmPosFlag = true;
      _rndmSphereCapRotMat = _ComputeRotMatrix(cos(_rndmSphereCapAlpha), sin(_rndmSphereCapAlpha),
                                               cos(_rndmSphereCapBeta), sin(_rndmSphereCapBeta));
      SetMinPhi(0);
      SetMaxPhi(360. * deg);
    }
    _ResetFlat();
  }
}

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

void GGSGunGeneratorAction::_SetRndmDirFlag() {
  if ((_rndmMinThetaFlag && _rndmMaxThetaFlag) || (_rndmMinPhiFlag && _rndmMaxPhiFlag))
    _rndmDirFlag = true;
  else
    _rndmDirFlag = false;
}

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

void GGSGunGeneratorAction::_SetRndmPosFlag() {
  if (_rndmMinPositionFlag && _rndmMaxPositionFlag)
    _rndmPosFlag = true;
  else
    _rndmPosFlag = false;
}

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

void GGSGunGeneratorAction::GeneratePrimaries(G4Event *anEvent) {

  static const std::string routineName("GGSGunGeneratorAction::GeneratePrimaries");
  static const G4ThreeVector xDirection(1, 0, 0);
  static const G4ThreeVector yDirection(0, 1, 0);

  _nDiscarded = 0;

  // Initialize shooting parameters with default values
  G4ThreeVector position(_position);
  G4ThreeVector direction(_direction);
  G4double energy(_energy);

  bool acceptanceCheck = true;
  do { // Acceptance check loop

    // ++++++++++++++ Randomize direction ++++++++++++++
    // The random direction computed here is always referred to the shooting frame.
    // Since the shooting frame may depend on shooting point (eg., shooting
    // from a sphere) the transformation from shooting frame to absolute
    // frame is done after random position have been chosen.
    if (_rndmDirFlag) {
      G4double cosPhi(_cosPhi);
      G4double sinPhi(_sinPhi);
      if (_rndmMinPhiFlag && _rndmMaxPhiFlag) {
        double phi = CLHEP::RandFlat::shoot(_minPhi, _maxPhi); // (0,2pi) or (-2pi,2pi) if shooting along -Z or -X
        cosPhi = cos(phi);
        sinPhi = sin(phi);
      }

      G4double sinTheta(_sinTheta);
      G4double cosTheta(_cosTheta);
      if (_rndmMinThetaFlag && _rndmMaxThetaFlag) {
        if (_rndmFromFlat || _rndmFromSphere) {
          // To generate an isotropic flux from a flat surface, the random angle must be uniform in cos^2 theta
          float cos2Theta = CLHEP::RandFlat::shoot(_minCos2ThetaRand, _maxCos2ThetaRand);
          cosTheta = sqrt(cos2Theta);
          sinTheta = sqrt(1 - cos2Theta);
        } else {
          // This will sample a uniform theta distribution on the sphere. This formula comes from the inverse CDF
          // method. See eg. http://en.wikibooks.org/wiki/Mathematica/Uniform_Spherical_Distribution
          // The form reported here is slightly different but equivalent (it can be easily seen using the
          // identity (sin x/2)^2 = (1 - cos x)/2 ).
          cosTheta = 1 - 2 * CLHEP::RandFlat::shoot(_minThetaRand, _maxThetaRand);
          sinTheta = sqrt(1 - cosTheta * cosTheta);
        }
      }

      // Compute position std::vector and direction versor (Z is shoot axis in the shooting frame)
      direction[0] = -sinTheta * cosPhi; // momentum direction (unitary vectors)
      direction[1] = -sinTheta * sinPhi;
      direction[2] = -cosTheta;
    }

    // ++++++++++++++ Randomize position ++++++++++++++
    if (_rndmPosFlag) {
      if (_rndmMinPositionFlag && _rndmMaxPositionFlag) {
        // This works for all "flat" domains (point, line, plane, parallelepiped)
        position[0] = CLHEP::RandFlat::shoot(_minPosition(0), _maxPosition(0));
        position[1] = CLHEP::RandFlat::shoot(_minPosition(1), _maxPosition(1));
        position[2] = CLHEP::RandFlat::shoot(_minPosition(2), _maxPosition(2));
        if (_rndmFromFlat) { // Transform shoot direction to absolute reference frame
          // Transform direction if shoot axis is not the Z axis.
          // 3rd coordinate is relative to shooting axis so it is x when shooting along X,
          // so the coordinate transformation to go to absolute reference frame is:
          // (y, z, x) -> (x, y, z) when shooting along X
          // (z, x, y) -> (x, y, z) when shooting along X
          if (_shootAxis == xDirection) {
            double buf = direction[0];
            direction[0] = direction[2];
            direction[2] = direction[1];
            direction[1] = buf;
          } else if (_shootAxis == yDirection) {
            double buf = direction[0];
            direction[0] = direction[1];
            direction[1] = direction[2];
            direction[2] = buf;
          }
        }
      } else if (_rndmFromSphere) {
        // Randomize position on the spherical cap (in the reference frame of the generation sphere).
        // Gamma and delta are the polar and azimuth angle of the position vector w.r.t. the
        // sphere center.
        float cosGamma = 1 - 2 * CLHEP::RandFlat::shoot(0., _rndmSphereCapExtRand);
        float delta = CLHEP::RandFlat::shoot(0., 2 * CLHEP::pi);
        float sinGamma = sqrt(1 - cosGamma * cosGamma);
        float sinDelta = sin(delta);
        float cosDelta = cos(delta);

        // Position vector in cap frame
        position[0] = _rndmSphereRadius * sinGamma * cosDelta;
        position[1] = _rndmSphereRadius * sinGamma * sinDelta;
        position[2] = _rndmSphereRadius * cosGamma;

        // Position vector in absolute (world) frame
        position *= _rndmSphereCapRotMat;
        position += _rndmSphereCenter;

        // Direction vector in absolute (world) frame
        // At this point, the direction vector is built so that the distribution of its polar angle with
        // respect to the normal incidence is flat in cos^2 (isotropic flux). The vector coordinates have
        // been constructed in the reference frame of the generation point on the sphere, so they must be
        // rotated to the reference frame of the cap center and then to the world reference frame, in this
        // order.
        direction *= _ComputeRotMatrix(cosGamma, sinGamma, cosDelta, sinDelta); // Direction
        direction *= _rndmSphereCapRotMat;
      }
    }

    // ++++++++++++++ Check acceptance ++++++++++++++
    if (_acceptanceFlag && (_rndmPosFlag || _rndmDirFlag)) {
      GGSDetectorConstruction *detConst =
          (GGSDetectorConstruction *)(GGSRunManager::GetRunManager()->GetUserDetectorConstruction());
      acceptanceCheck = detConst->GetGeometry()->IsInsideAcceptance(position, direction);
      if (!acceptanceCheck)
        _nDiscarded++;
    }

  } while (!acceptanceCheck);

  // ++++++++++++++ Randomize energy ++++++++++++++
  if (_rndmMinEnergyFlag && _rndmMaxEnergyFlag) {
    if (_spectralIndexFlag) {
      energy = _GenSpectrumPowerLaw(_minEnergy, _maxEnergy, _spectralIndex);
    } else {
      energy = CLHEP::RandFlat::shoot(_minEnergy, _maxEnergy);
    }
  }

  COUT(DEEPDEB) << "Generated primary:" << ENDL;
  CCOUT(DEEPDEB) << "  PDG code: " << _gunGenerator->GetParticleDefinition()->GetPDGEncoding() << ENDL;
  CCOUT(DEEPDEB) << "  position: " << position << ENDL;
  CCOUT(DEEPDEB) << "  direction: " << direction << ENDL;
  CCOUT(DEEPDEB) << "  energy: " << energy << ENDL;

  //
  // Set values for primary particle and shoot
  //
  _gunGenerator->SetParticlePosition(position);
  _gunGenerator->SetParticleMomentumDirection(direction);
  _gunGenerator->SetParticleEnergy(energy);

  // Shoot
  _gunGenerator->GeneratePrimaryVertex(anEvent);
}

GGSParameters GGSGunGeneratorAction::GetParameters() const {
  GGSParameters params;
  params.SetParam("generator", std::string("gun"));

  params.SetParam("acceptanceCheck", _acceptanceFlag);

  // Particle
  params.SetParam("particlePdgCode", _gunGenerator->GetParticleDefinition()->GetPDGEncoding());

  // Position
  if (_rndmPosFlag) {
    if (_rndmMinPositionFlag && _rndmMaxPositionFlag) {
      params.SetParam("minPosition", std::vector<double>{_minPosition[0] / CLHEP::cm, _minPosition[1] / CLHEP::cm,
                                                         _minPosition[2] / CLHEP::cm});
      params.SetParam("maxPosition", std::vector<double>{_maxPosition[0] / CLHEP::cm, _maxPosition[1] / CLHEP::cm,
                                                         _maxPosition[2] / CLHEP::cm});
      if (_rndmFromFlat) {
        params.SetParam("genPosition", std::string("plane"));
      } else {
        params.SetParam("genPosition", std::string("box"));
      }
    } else if (_rndmFromSphere) {
      params.SetParam("genPosition", std::string("sphere"));
      params.SetParam("sphereCenter",
                      std::vector<double>{_rndmSphereCenter[0] / CLHEP::cm, _rndmSphereCenter[1] / CLHEP::cm,
                                          _rndmSphereCenter[2] / CLHEP::cm});
      params.SetParam("sphereCapPosition",
                      std::vector<double>{_rndmSphereCapPosition[0] / CLHEP::cm, _rndmSphereCapPosition[1] / CLHEP::cm,
                                          _rndmSphereCapPosition[2] / CLHEP::cm});
      params.SetParam("sphereCapExtension", acos(-2. * _rndmSphereCapExtRand + 1) / CLHEP::rad);
    }
  } else {
    params.SetParam("genPosition", std::string("point"));
    params.SetParam("position",
                    std::vector<double>{_position[0] / CLHEP::cm, _position[1] / CLHEP::cm, _position[2] / CLHEP::cm});
  }

  // Direction
  if (_rndmDirFlag) {
    bool isotropic = false;
    if (_rndmMinPhiFlag && _rndmMaxPhiFlag) {
      params.SetParam("minPhi", _minPhi / CLHEP::rad);
      params.SetParam("maxPhi", _maxPhi / CLHEP::rad);
    }
    if (_rndmMinThetaFlag && _rndmMaxThetaFlag) {
      params.SetParam("minTheta", acos(-2. * _minThetaRand + 1) / CLHEP::rad);
      params.SetParam("maxTheta", acos(-2. * _maxThetaRand + 1) / CLHEP::rad);
      if (_rndmFromFlat || _rndmFromSphere) {
        isotropic = true;
      }
    }
    if (isotropic) {
      params.SetParam("genDirection", std::string("isotropic"));
    } else {
      params.SetParam("genDirection", std::string("uniform"));
    }
  } else {
    params.SetParam("genDirection", std::string("fixed"));
    params.SetParam("theta", acos(_cosTheta) / CLHEP::rad);
    params.SetParam("phi", acos(_cosPhi) / CLHEP::rad);
  }

  // Energy
  if (_rndmMinEnergyFlag && _rndmMaxEnergyFlag) {
    params.SetParam("minEnergy", _minEnergy / CLHEP::GeV);
    params.SetParam("maxEnergy", _maxEnergy / CLHEP::GeV);
    if (_spectralIndexFlag) {
      params.SetParam("genSpectrum", std::string("powerLaw"));
      params.SetParam("spectralIndex", _spectralIndex);
    } else {
      params.SetParam("genSpectrum", std::string("flat"));
    }
  } else {
    params.SetParam("genSpectrum", std::string("monochromatic"));
    params.SetParam("energy", _energy / CLHEP::GeV);
  }

  return params;
}

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

// Generate spectrum power law f(E)= E^gamma (when gamma==0 then is uniform generation).
G4double GGSGunGeneratorAction::_GenSpectrumPowerLaw(G4double Emin, G4double Emax, G4double gamma) {
  G4double energy;
  // Uniform generation
  if (gamma == 0.) {
    energy = CLHEP::RandFlat::shoot(Emin, Emax);
    return energy;
  }
  G4double alpha = 1. + gamma; // integral spectral index
  if (alpha == 0.) {
    energy = exp(log(Emin) + CLHEP::RandFlat::shoot(0., 1.) * (log(Emax) - log(Emin)));
  } else {
    if (Emin == 0.)
      Emin = 1.E-10;
    energy =
        pow((CLHEP::RandFlat::shoot(0., 1.) * (pow(Emax, alpha) - pow(Emin, alpha)) + pow(Emin, alpha)), 1. / alpha);
  }

  return energy;
}