Commit 829edf0e authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Merge branch 'DepositionGeant4_RectangleBeam' into 'master'

DepositionGeant4: Allow for a beam of rectangular and elliptical shape.

See merge request allpix-squared/allpix-squared!1077
parents 5cb9cc94 3a33e3a0
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -26,6 +26,7 @@ The following authors, in alphabetical order, have developed or contributed to A
* Carsten Daniel Burgard, DESY, [cburgard](https://gitlab.cern.ch/cburgard)
* Maximilian Felix Caspar, DESY, [mcaspar](https://gitlab.cern.ch/mcaspar)
* Liejian Chen, Institute of High Energy Physics Beijing, [chenlj](https://github.com/chenlj)
* Naomi Davis, DESY, [naomi](https://gitlab.cern.ch/naomi)
* Manuel Alejandro Del Rio Viera, DESY, [mdelriov](https://gitlab.cern.ch/mdelriov)
* Katharina Dort, University of Gießen, [kdort](https://gitlab.cern.ch/kdort)
* Neal Gauvin, Université de Genève, [ngauvin](https://gitlab.cern.ch/ngauvin)
+25 −0
Original line number Diff line number Diff line
@@ -407,6 +407,10 @@ void DepositionGeant4Module::run(Event* event) {

                double deposited_energy = static_cast<double>(Units::convert(sensor->getDepositedEnergy(), "keV"));
                energy_per_event_[sensor->getName()]->Fill(deposited_energy);

                for(const auto& track_position : sensor->getTrackIncidentPositions()) {
                    incident_track_position_[sensor->getName()]->Fill(track_position.x(), track_position.y());
                }
            }
        }
    } catch(AbortEventException& e) {
@@ -432,6 +436,9 @@ void DepositionGeant4Module::finalize() {
        for(auto& histogram : energy_per_event_) {
            histogram.second->Write();
        }
        for(auto& histogram : incident_track_position_) {
            histogram.second->Write();
        }
    }

    // Print summary or warns if module did not output any charges
@@ -536,6 +543,10 @@ void DepositionGeant4Module::construct_sensitive_detectors_and_fields() {
                (static_cast<int>(maximum_charge / 2. * Units::convert(charge_creation_energy, "eV")) / 10) * 10 + 10;
            int nbins = 5 * maximum_charge;

            // Get detector model size
            auto sensor_size = detector->getModel()->getSensorSize();
            auto pixel_size = detector->getModel()->getPixelSize();

            // Create histograms if needed
            {
                std::lock_guard<std::mutex> lock(histogram_mutex_);
@@ -560,6 +571,20 @@ void DepositionGeant4Module::construct_sensitive_detectors_and_fields() {
                                              0,
                                              maximum_energy);
                }

                plot_name = "incident_track_position_" + sensitive_detector_action->getName();

                if(incident_track_position_.find(sensitive_detector_action->getName()) == incident_track_position_.end()) {
                    incident_track_position_[sensitive_detector_action->getName()] =
                        CreateHistogram<TH2D>(plot_name.c_str(),
                                              "incident track position;X [mm];Y [mm];Z",
                                              5000,
                                              -pixel_size.X() / 2,
                                              sensor_size.X() - pixel_size.X() / 2,
                                              5000,
                                              -pixel_size.Y() / 2,
                                              sensor_size.Y() - pixel_size.Y() / 2);
                }
            }
        }
    }
+2 −0
Original line number Diff line number Diff line
@@ -30,6 +30,7 @@
#include "tools/ROOT.h"

#include <TH1D.h>
#include <TH2D.h>

class G4UserLimits;
class G4RunManager;
@@ -122,6 +123,7 @@ namespace allpix {
        // Vector of histogram pointers for debugging plots
        std::map<std::string, Histogram<TH1D>> charge_per_event_;
        std::map<std::string, Histogram<TH1D>> energy_per_event_;
        std::map<std::string, Histogram<TH2D>> incident_track_position_;

        // Total deposited charges
        std::atomic_uint total_charges_{0};
+39 −5
Original line number Diff line number Diff line
@@ -22,10 +22,12 @@
#include <G4ParticleTable.hh>
#include <G4RunManager.hh>
#include <G4UImanager.hh>
#include <Math/Vector2D.h>
#include <core/module/exceptions.h>

#include "core/config/exceptions.h"
#include "core/utils/log.h"
#include "tools/ROOT.h"
#include "tools/geant4/geant4.h"

using namespace allpix;
@@ -127,13 +129,45 @@ GeneratorActionG4::GeneratorActionG4(const Configuration& config)

            // Set position parameters
            single_source->GetPosDist()->SetPosDisType("Beam");
            auto beam_size = config.get<double>("beam_size", 0);
            if(config.get<bool>("flat_beam", false) == true) {

            // Get beam_size parameter(s) from config file
            ROOT::Math::XYVector beam_size{};
            try {
                beam_size = config_.get<ROOT::Math::XYVector>("beam_size", {0, 0});
            } catch(InvalidKeyError&) {
                const auto size = config_.get<double>("beam_size", 0);
                beam_size = {size, size};
            }

            auto beam_shape = config_.get<BeamShape>("beam_shape", BeamShape::CIRCLE);
            if(config_.get<bool>("flat_beam", false)) {
                // Note: G4 definition of rectangle swaps x and y wrt our global coordinate system
                single_source->GetPosDist()->SetPosDisType("Plane");
                if(beam_shape == BeamShape::RECTANGLE) {
                    single_source->GetPosDist()->SetPosDisShape("Rectangle");
                    single_source->GetPosDist()->SetHalfX((beam_size.y()) / 2);
                    single_source->GetPosDist()->SetHalfY((beam_size.x()) / 2);
                } else if(beam_shape == BeamShape::CIRCLE) {
                    single_source->GetPosDist()->SetPosDisShape("Circle");
                single_source->GetPosDist()->SetRadius(beam_size);
                    single_source->GetPosDist()->SetRadius(beam_size.x());
                } else if(beam_shape == BeamShape::ELLIPSE) {
                    single_source->GetPosDist()->SetPosDisShape("Ellipse");
                    single_source->GetPosDist()->SetHalfX((beam_size.y()) / 2);
                    single_source->GetPosDist()->SetHalfY((beam_size.x()) / 2);
                } else {
                single_source->GetPosDist()->SetBeamSigmaInR(beam_size);
                    throw InvalidValueError(config_, "beam_shape", "Cannot use this beam shape with flat beams");
                }
            } else {
                if(beam_shape == BeamShape::CIRCLE) {
                    single_source->GetPosDist()->SetBeamSigmaInR(beam_size.x());
                } else if(beam_shape == BeamShape::ELLIPSE || beam_shape == BeamShape::RECTANGLE) {
                    single_source->GetPosDist()->SetBeamSigmaInX((beam_size.y()) / 2);
                    single_source->GetPosDist()->SetBeamSigmaInY((beam_size.x()) / 2);
                } else {
                    throw InvalidValueError(config_, "beam_shape", "This beam shape can only be used with flat beams");
                }
            }

            single_source->GetPosDist()->SetPosRot1(angref1);
            single_source->GetPosDist()->SetPosRot2(angref2);

+9 −0
Original line number Diff line number Diff line
@@ -40,6 +40,15 @@ namespace allpix {
            POINT,  ///< Point source
        };

        /**
         * @brief Different shapes of particle beams
         */
        enum class BeamShape {
            CIRCLE,    ///< Circular beam
            ELLIPSE,   ///< Elliptical beam
            RECTANGLE, ///< Rectangular beam
        };

        /**
         * @brief Constructs the generator action
         * @param config Configuration of the \ref DepositionGeant4Module module
Loading