ghit.h

Go to the documentation of this file.

#pragma once
 
// geant4
#include "G4VHit.hh"
#include "G4THitsCollection.hh"
#include "G4Allocator.hh"
 
#include "G4ThreeVector.hh"
#include "G4Step.hh"
#include "G4Colour.hh"
 
// gemc
#include <gemc/gtouchable/gtouchable.h>
 
// c++
#include <optional>
#include <atomic>
#include <map>
 
class GHit : public G4VHit
{
public:
    GHit(std::shared_ptr<GTouchable> gt, const G4Step* thisStep = nullptr,
         const std::string&          cScheme                    = "default");
 
    ~GHit() override = default;
 
    inline void* operator new(size_t);
 
    inline void operator delete(void*);
 
    void Draw() override;
 
    [[nodiscard]] bool is_same_hit(const GHit* hit) const;
 
private:
    G4Colour colour_touch, colour_hit, colour_passby;
 
    bool setColorSchema();
 
    std::string colorSchema;
 
    std::shared_ptr<GTouchable> gtouchable;
 
    // -------------------------------------------------------------------------
    // Per-step data (vectors)
    // -------------------------------------------------------------------------
 
    std::vector<double> edeps;
 
    std::vector<double> times;
 
    std::vector<G4ThreeVector> globalPositions;
 
    std::vector<G4ThreeVector> localPositions;
 
    std::vector<G4ThreeVector> trackVertexPositions;
 
    std::vector<G4ThreeVector> motherTrackVertexPositions;
 
    std::vector<int> pids;
 
    std::vector<int> tids;
 
    std::vector<int> motherTids;
 
    std::vector<std::string> processNames;
 
    std::vector<G4ThreeVector> momenta;
 
    std::vector<double> trackEs;
 
    std::vector<int> motherPids;
 
    static thread_local std::map<int, int> pdgById;
 
    // -------------------------------------------------------------------------
    // Aggregated / calculated quantities (lazy)
    // -------------------------------------------------------------------------
 
    std::optional<double> totalEnergyDeposited;
 
 
    double averageTime;
 
    G4ThreeVector avgGlobalPosition;
 
    G4ThreeVector avgLocalPosition;
 
    std::string processName;
 
 
    static std::atomic<int> globalHitCounter;
 
    static thread_local std::map<int, G4ThreeVector> trackVertexById;
 
public:
    // -------------------------------------------------------------------------
    // Inline accessors (returning copies by design)
    // -------------------------------------------------------------------------
 
    [[nodiscard]] inline std::vector<double> getEdeps() const { return edeps; }
 
    [[nodiscard]] inline std::vector<double> getTimes() const { return times; }
 
    [[nodiscard]] inline std::vector<G4ThreeVector> getGlobalPositions() const { return globalPositions; }
 
    [[nodiscard]] inline std::vector<G4ThreeVector> getLocalPositions() const { return localPositions; }
 
    [[nodiscard]] inline std::vector<G4ThreeVector> getTrackVertexPositions() const { return trackVertexPositions; }
 
    [[nodiscard]] inline G4ThreeVector getTrackVertexPosition() const { return trackVertexPositions.front(); }
 
    [[nodiscard]] inline std::vector<G4ThreeVector> getMotherTrackVertexPositions() const {
        return motherTrackVertexPositions;
    }
 
    [[nodiscard]] inline G4ThreeVector getMotherTrackVertexPosition() const {
        return motherTrackVertexPositions.front();
    }
 
    [[nodiscard]] inline std::vector<int> getPids() const { return pids; }
 
    [[nodiscard]] inline int getPid() const { return pids.front(); }
 
    [[nodiscard]] inline  std::vector<int> getTids() const { return tids; }
 
    [[nodiscard]] inline int getTid() const { return tids.front(); }
 
    [[nodiscard]] inline std::vector<int> getMotherTids() const { return motherTids; }
 
    [[nodiscard]] inline int getMotherTid() const { return motherTids.front(); }
 
    [[nodiscard]] inline double getE() const { return trackEs.front(); }
 
    [[nodiscard]] inline size_t nsteps() const { return edeps.size(); }
 
    [[nodiscard]] inline size_t getStepCount() const { return edeps.size(); }
 
    [[nodiscard]] inline std::vector<G4ThreeVector> getMomenta() const { return momenta; }
 
    [[nodiscard]] inline G4ThreeVector getMomentum() const { return momenta.front(); }
 
    [[nodiscard]] inline std::vector<double> getTrackEs() const { return trackEs; }
 
    [[nodiscard]] inline double getTrackE() const { return trackEs.front(); }
 
    [[nodiscard]] inline std::vector<int> getMotherPids() const { return motherPids; }
 
    [[nodiscard]] inline int getMpid() const { return motherPids.front(); }
 
    [[nodiscard]] inline std::string getProcID() const { return processName; }
 
    [[nodiscard]] inline std::string getProcessName() const { return processName; }
 
    [[nodiscard]] inline std::shared_ptr<GTouchable> getGTouchable() const { return gtouchable; }
 
    [[nodiscard]] inline std::vector<GIdentifier> getGID() const { return gtouchable->getIdentity(); }
 
    [[nodiscard]] inline std::vector<double> getDetectorDimensions() const {
        return gtouchable->getDetectorDimensions();
    }
 
    [[nodiscard]] inline double getMass() const { return gtouchable->getMass(); }
 
    // -------------------------------------------------------------------------
    // Aggregation / calculation API
    // -------------------------------------------------------------------------
 
    void calculateInfos();
 
    double getTotalEnergyDeposited();
 
    double getAverageTime();
 
    G4ThreeVector getAvgLocalPosition();
 
    G4ThreeVector getAvgGlobaPosition();
 
    // -------------------------------------------------------------------------
    // Hit filling / testing helpers
    // -------------------------------------------------------------------------
 
    void addHitInfos(const G4Step* thisStep);
 
    static void clearTrackVertexCache();
 
    void randomizeHitForTesting(int nsteps);
 
    [[nodiscard]] std::vector<int> getTTID() const;
 
    static GHit* create(const std::shared_ptr<GOptions>& gopts) {
        auto gt  = GTouchable::create(gopts);
        auto hit = new GHit(gt);
        // Randomize between 1 and 10 steps in a deterministic, thread-safe manner.
        hit->randomizeHitForTesting(1 + globalHitCounter.fetch_add(1, std::memory_order_relaxed) % 10);
        return hit;
    }
};
 
// MT definitions, as from:
// https://twiki.cern.ch/twiki/bin/view/Geant4/QuickMigrationGuideForGeant4V10
extern G4ThreadLocal G4Allocator<GHit>* GHitAllocator;
using GHitsCollection = G4THitsCollection<GHit>;
 
inline void* GHit::operator new(size_t) {
    if (!GHitAllocator) GHitAllocator = new G4Allocator<GHit>;
    return (void*)GHitAllocator->MallocSingle();
}
 
inline void GHit::operator delete(void* hit) {
    if (!GHitAllocator) { GHitAllocator = new G4Allocator<GHit>; }
 
    GHitAllocator->FreeSingle((GHit*)hit);
}
 
[[nodiscard]] inline std::string getIdentityString(std::vector<GIdentifier> gidentity) {
    // Build a compact label from the stored identifier vector.
    std::string identifierString;
    for (size_t i = 0; i < gidentity.size() - 1; i++) {
        identifierString += gidentity[i].getName() + "->" + std::to_string(gidentity[i].getValue()) + ", ";
    }
    identifierString += gidentity.back().getName() + "->" + std::to_string(gidentity.back().getValue());
    return identifierString;
}
 
[[nodiscard]] inline std::map<std::string, int> getIdentityMap(std::vector<GIdentifier> gidentity) {
    std::map<std::string, int> identityMap;
    for (auto& id : gidentity) {
        identityMap[id.getName()] = id.getValue();
    }
    return identityMap;
}