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:
    struct MotherInfo
    {
        int                          trackId{};
        std::optional<G4ThreeVector> vertex;
        std::optional<int>           pid;
    };

    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<MotherInfo> motherInfos;

    std::vector<int> pids;

    std::vector<int> tids;

    std::vector<std::string> processNames;

    std::vector<G4ThreeVector> momenta;

    std::vector<double> trackEs;

    static thread_local std::map<int, int> pdgById;
 
    // -------------------------------------------------------------------------
    // Aggregated / calculated quantities (lazy)
    // -------------------------------------------------------------------------

    struct CalculatedState {
        double                     totalEnergyDeposited{};
        double                     averageTime{};
        G4ThreeVector              averageGlobalPosition;
        G4ThreeVector              averageLocalPosition;
        std::optional<std::string> processName;
    };

    mutable std::optional<CalculatedState> calculatedState;

    [[nodiscard]] const CalculatedState& getCalculatedState() const;

    void invalidateCalculatedState();
 

    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 const std::vector<MotherInfo>& getMotherInfos() const { return motherInfos; }

    [[nodiscard]] inline const MotherInfo& getMotherInfo() const { return motherInfos.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 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]] size_t getNumberOfOpticalPhotons() const;

    [[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]] const std::optional<std::string>& getProcessName() const;

    [[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
    // -------------------------------------------------------------------------

    [[nodiscard]] double getTotalEnergyDeposited() const;

    [[nodiscard]] double getAverageTime() const;

    [[nodiscard]] G4ThreeVector getAvgLocalPosition() const;

    [[nodiscard]] G4ThreeVector getAvgGlobalPosition() const;
 
    // -------------------------------------------------------------------------
    // 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;
}