Commit 121c6e8a authored by Naomi Davis's avatar Naomi Davis Committed by Simon Spannagel
Browse files

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

parent 4fa539be
Loading
Loading
Loading
Loading
+21 −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(auto track_position : sensor->getIncidentPosition()){
                    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
@@ -560,6 +567,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",
                                              2000, 
                                              -50, 
                                              50, 
                                              6000, 
                                              -30, 
                                              30);
                }
            }
        }
    }
+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};
+32 −5
Original line number Diff line number Diff line
@@ -23,9 +23,16 @@
#include <G4RunManager.hh>
#include <G4UImanager.hh>
#include <core/module/exceptions.h>
#include <Math/Point2D.h>
#include <Math/Vector2D.h>
#include <Math/Translation3D.h>


#include "TMath.h"

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

using namespace allpix;
@@ -127,12 +134,32 @@ 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");   
            } catch(InvalidKeyError &) {
                auto vector_size = config_.get<double>("beam_size", 0);   
                beam_size = {vector_size, vector_size};  
            }

            if(config_.get<bool>("flat_beam", false) == true) {
                if(config_.get<std::string>("beam_shape", "Circle") == "Rectangle") {
                    single_source->GetPosDist()->SetPosDisType("Plane"); //?
                    single_source->GetPosDist()->SetPosDisShape("Rectangle");
                    single_source->GetPosDist()->SetHalfX((beam_size.x()) / 2);
                    single_source->GetPosDist()->SetHalfY((beam_size.y()) / 2);
                } else {
                    single_source->GetPosDist()->SetPosDisShape("Circle");
                single_source->GetPosDist()->SetRadius(beam_size);
                    single_source->GetPosDist()->SetRadius(beam_size.x());
                }
            } else {
                single_source->GetPosDist()->SetBeamSigmaInR(beam_size);
                if(config_.get<std::string>("beam_shape", "Circle") == "Ellipse"){
                    single_source->GetPosDist()->SetBeamSigmaInX((beam_size.x()) / 2);
                    single_source->GetPosDist()->SetBeamSigmaInY((beam_size.y()) / 2); 
                }
                single_source->GetPosDist()->SetBeamSigmaInR(beam_size.x());
            }
            single_source->GetPosDist()->SetPosRot1(angref1);
            single_source->GetPosDist()->SetPosRot2(angref2);
+6 −0
Original line number Diff line number Diff line
@@ -136,6 +136,10 @@ double SensitiveDetectorActionG4::getTotalDepositedEnergy() const { return total

double SensitiveDetectorActionG4::getDepositedEnergy() const { return deposited_energy_; }

std::vector<ROOT::Math::XYZPoint> SensitiveDetectorActionG4::getIncidentPosition() const {    
    return incident_track_position_;
}

void SensitiveDetectorActionG4::clearEventInfo() {
    LOG(DEBUG) << "Clearing track and deposit vectors";

@@ -188,6 +192,8 @@ void SensitiveDetectorActionG4::dispatchMessages(Module* module, Messenger* mess
        mc_particles.back().setKineticEnergyStart(track_kinetic_energy_start_.at(track_id));
        id_to_particle_[track_id] = mc_particles.size() - 1;

        incident_track_position_.push_back(local_begin);

        LOG(DEBUG) << "Found MC particle " << pdg_code << " crossing detector " << detector_->getName() << " from "
                   << Units::display(local_begin, {"mm", "um"}) << " to " << Units::display(local_end, {"mm", "um"})
                   << " local after " << Units::display(track_time_global, {"ns", "ps"}) << " global / "
+7 −0
Original line number Diff line number Diff line
@@ -72,6 +72,11 @@ namespace allpix {
         */
        double getDepositedEnergy() const;

        /**
         * @brief Get the position of the incident particle tracks for this event.
         */
        std::vector<ROOT::Math::XYZPoint> getIncidentPosition() const;

        /**
         * @brief Clears depopsition information vectors in preparation for the next event.
         */
@@ -126,6 +131,8 @@ namespace allpix {

        // List of positions for deposits
        std::vector<ROOT::Math::XYZPoint> deposit_position_;
        std::vector<ROOT::Math::XYZPoint> incident_track_position_;

        std::vector<unsigned int> deposit_charge_;
        std::vector<double> deposit_energy_;
        std::vector<double> deposit_time_;