diff --git a/source/generators/LED450LWcone.cc b/source/generators/LED450LWcone.cc new file mode 100644 index 0000000000..60c72911fa --- /dev/null +++ b/source/generators/LED450LWcone.cc @@ -0,0 +1,160 @@ +// ---------------------------------------------------------------------------- +// nexus | LED450LWcone.cc +// +// This class is the primary generator for events consisting of +// a single particle. The user must specify via configuration +// parameters the particle type, a kinetic energy interval and, optionally, +// a momentum direction. +// Particle energy is generated with flat random probability +// between E_min and E_max. +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#include "LED450LWcone.h" + +#include "DetectorConstruction.h" +#include "GeometryBase.h" +#include "RandomUtils.h" +#include "FactoryBase.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "CLHEP/Units/SystemOfUnits.h" + +using namespace nexus; +using namespace CLHEP; + +REGISTER_CLASS(LED450LWcone, G4VPrimaryGenerator) + + +LED450LWcone::LED450LWcone(): +G4VPrimaryGenerator(), msg_(0), particle_definition_(0), +energy_min_(0.), energy_max_(0.), geom_(0), momentum_{}, +costheta_min_(-1.), costheta_max_(1.), phi_min_(0.), phi_max_(2.*pi) +{ + msg_ = new G4GenericMessenger(this, "/Generator/LED450LWcone/", + "Control commands of single-particle generator."); + + msg_->DeclareMethod("particle", &LED450LWcone::SetParticleDefinition, + "Set particle to be generated."); + + G4GenericMessenger::Command& min_energy = + msg_->DeclareProperty("min_energy", energy_min_, + "Set minimum kinetic energy of the particle."); + min_energy.SetUnitCategory("Energy"); + min_energy.SetParameterName("min_energy", false); + min_energy.SetRange("min_energy>0."); + + G4GenericMessenger::Command& max_energy = + msg_->DeclareProperty("max_energy", energy_max_, + "Set maximum kinetic energy of the particle"); + max_energy.SetUnitCategory("Energy"); + max_energy.SetParameterName("max_energy", false); + max_energy.SetRange("max_energy>0."); + + msg_->DeclareProperty("region", region_, + "Set the region of the geometry where the vertex will be generated."); + + + msg_->DeclarePropertyWithUnit("momentum", "mm", momentum_, "Set particle 3-momentum."); + + msg_->DeclareProperty("min_costheta", costheta_min_, + "Set minimum cosTheta for the direction of the particle."); + msg_->DeclareProperty("max_costheta", costheta_max_, + "Set maximum cosTheta for the direction of the particle."); + msg_->DeclareProperty("min_phi", phi_min_, + "Set minimum phi for the direction of the particle."); + msg_->DeclareProperty("max_phi", phi_max_, + "Set maximum phi for the direction of the particle."); + + + DetectorConstruction* detconst = (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + geom_ = detconst->GetGeometry(); +} + + + +LED450LWcone::~LED450LWcone() +{ + delete msg_; +} + + + +void LED450LWcone::SetParticleDefinition(G4String particle_name) +{ + particle_definition_ = + G4ParticleTable::GetParticleTable()->FindParticle(particle_name); + + if (!particle_definition_) + G4Exception("[LED450LW]", "SetParticleDefinition()", + FatalException, "User gave an unknown particle name."); +} + + + +void LED450LWcone::GeneratePrimaryVertex(G4Event* event) +{ + // Generate uniform random energy in [E_min, E_max] + G4double kinetic_energy = nexus::UniformRandomInRange(energy_max_, energy_min_); + + // Calculate cartesian components of momentum + G4double mass = particle_definition_->GetPDGMass(); + G4double energy = kinetic_energy + mass; + G4double pmod = std::sqrt(energy*energy - mass*mass); + + bool fixed_momentum = momentum_ != G4ThreeVector{}; + bool restrict_angle = costheta_min_ != -1. || costheta_max_ != 1. || phi_min_ != 0. || phi_max_ !=2.*pi; + + G4ThreeVector dir; // it will be set in the if branches below + if (fixed_momentum) { // if the user provides a momentum direction + dir = momentum_.unit(); + } else if (restrict_angle) { // if the user provides a range of angles + dir = RandomDirectionInRange(costheta_min_, costheta_max_, phi_min_, phi_max_); + } + + G4ThreeVector p_dir = G4LambertianRand(dir); + + G4double coseno = -1*p_dir.dot(momentum_); + + G4ThreeVector p; + + while((costheta_min_ > coseno) || (coseno > costheta_max_)){ + p_dir = G4LambertianRand(dir); + coseno = -1*p_dir.dot(dir); + } + + p = pmod * p_dir; + + // Create the new primary particle and set it some properties + auto particle = new G4PrimaryParticle(particle_definition_, p.x(), p.y(), p.z()); + + // Set random polarization + if (particle_definition_ == G4OpticalPhoton::Definition()) { + G4ThreeVector polarization = G4RandomDirection(); + particle->SetPolarization(polarization); + } + + // Generate an initial position for the particle using the geometry + G4ThreeVector position = geom_->GenerateVertex(region_); + + // Particle generated at start-of-event + G4double time = 0.; + + // Create a new vertex + G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time); + + // Add particle to the vertex and this to the event + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); +} diff --git a/source/generators/LED450LWcone.h b/source/generators/LED450LWcone.h new file mode 100644 index 0000000000..4a05fda9d5 --- /dev/null +++ b/source/generators/LED450LWcone.h @@ -0,0 +1,73 @@ +// ---------------------------------------------------------------------------- +// nexus | LED450LWcone.h +// +// This class is the primary generator for events consisting of +// a single particle. The user must specify via configuration +// parameters the particle type, a kinetic energy interval and, optionally, +// a momentum direction. +// Particle energy is generated with flat random probability +// between E_min and E_max. +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#ifndef LED450LW_CONE_H +#define LED450LW_CONE_H + +#include + +class G4GenericMessenger; +class G4Event; +class G4ParticleDefinition; + + +namespace nexus { + + class GeometryBase; + + class LED450LWcone: public G4VPrimaryGenerator + { + public: + /// Constructor + LED450LWcone(); + /// Destructor + ~LED450LWcone(); + + /// This method is invoked at the beginning of the event. It sets + /// a primary vertex (that is, a particle in a given position and time) + /// in the event. + void GeneratePrimaryVertex(G4Event*); + + private: + + void SetParticleDefinition(G4String); + + /// Generate a random kinetic energy with flat probability in + // the interval [energy_min, energy_max]. + G4double RandomEnergy() const; + + private: + G4GenericMessenger* msg_; + + G4ParticleDefinition* particle_definition_; + + G4double energy_min_; ///< Minimum kinetic energy + G4double energy_max_; ///< Maximum kinetic energy + + const GeometryBase* geom_; ///< Pointer to the detector geometry + + G4String region_; + + G4ThreeVector momentum_; + + G4double costheta_min_; + G4double costheta_max_; + G4double phi_min_; + G4double phi_max_; + + + }; + +} // end namespace nexus + +#endif