geant4用G4Scintillation获取粒子打在闪烁体上的波形，并将输出的波形数据保存在txt文件中，用代码演示一下

以下是一个简单的示例代码，演示如何使用Geant4中的G4Scintillation模型获取闪烁体上的光信号，并将波形数据保存在txt文件中：

#include "G4Scintillation.hh"
#include "G4Scintillation.hh"
#include "G4OpticalPhoton.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
#include "G4EmProcessSubType.hh"
#include "G4OpBoundaryProcess.hh"
#include "G4RunManager.hh"
#include "G4Event.hh"
#include "G4SDManager.hh"
#include "G4HCofThisEvent.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include "G4ThreeVector.hh"
#include "G4ios.hh"
#include <fstream>

int main(int argc, char** argv)
{
    // Initialize Geant4 run manager
    G4RunManager* runManager = new G4RunManager;

    // Create a scintillation process
    G4Scintillation* scintProcess = new G4Scintillation();
    scintProcess->SetScintillationYieldFactor(1.0);
    scintProcess->SetTrackSecondariesFirst(true);
    scintProcess->SetVerboseLevel(0);

    // Register the scintillation process to the particle table
    G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
    G4ParticleDefinition* particle = particleTable->FindParticle("e-");
    G4EmProcessSubType processType = fScintillation;
    G4ProcessManager* pmanager = particle->GetProcessManager();
    pmanager->AddProcess(scintProcess, -1, 1, 1, processType);

    // Initialize the run manager
    runManager->Initialize();

    // Create an optical photon
    G4ParticleDefinition* opticalPhoton = G4OpticalPhoton::OpticalPhotonDefinition();

    // Create an event and a primary particle (an electron)
    G4Event* event = new G4Event();
    G4PrimaryVertex* vertex = new G4PrimaryVertex(G4ThreeVector(0, 0, 0), 0);
    G4PrimaryParticle* particle = new G4PrimaryParticle(G4Electron::ElectronDefinition(), 1.0, 0.0, 0.0, G4ThreeVector(0, 0, -1));
    vertex->SetPrimary(particle);
    event->AddPrimaryVertex(vertex);

    // Create a sensitive detector for the scintillation photons
    G4SDManager* sdManager = G4SDManager::GetSDMpointer();
    G4String sdName = "ScintillatorSD";
    ScintillatorSD* scintSD = new ScintillatorSD(sdName);
    sdManager->AddNewDetector(scintSD);

    // Set up the optical boundary processes
    G4OpBoundaryProcess* opBoundaryProcess = new G4OpBoundaryProcess();
    opBoundaryProcess->SetVerboseLevel(0);

    // Generate the scintillation photons
    runManager->SetUserAction(new GenScintillationPhotons(event, opticalPhoton, scintProcess, opBoundaryProcess, scintSD));

    // Run the simulation
    runManager->BeamOn(1);

    // Save the waveforms to a txt file
    std::ofstream outfile("waveform.txt");
    for (int i = 0; i < scintSD->GetWaveformSize(); i++) {
        outfile << scintSD->GetWaveformValue(i) << std::endl;
    }
    outfile.close();

    // Clean up
    delete runManager;
    delete event;
    delete vertex;
    delete particle;
    delete opticalPhoton;
    delete scintProcess;
    delete opBoundaryProcess;
    delete scintSD;

    return 0;
}

其中，`GenScintillationPhotons` 类是一个用户定义的类，用于在Geant4模拟中生成闪烁体上的光子，并将它们传递给一个敏感探测器。下面是这个类的代码：

class GenScintillationPhotons : public G4UserEventAction
{
public:
    GenScintillationPhotons(G4Event* anEvent, G4ParticleDefinition* optPhotonDef, G4Scintillation* scintProcess, G4OpBoundaryProcess* opBoundaryProcess, ScintillatorSD* scintSD)
        : G4UserEventAction(), fEvent(anEvent), fOpticalPhoton(optPhotonDef), fScintProcess(scintProcess), fOpBoundaryProcess(opBoundaryProcess), fScintSD(scintSD)
    {}

    virtual ~GenScintillationPhotons() {}

    virtual void BeginOfEventAction(const G4Event* anEvent)
    {
        fScintSD->InitializeWaveform();
    }

    virtual void EndOfEventAction(const G4Event* anEvent)
    {
        // Get the list of tracks in the event
        G4TrajectoryContainer* trajectoryContainer = anEvent->GetTrajectoryContainer();
        std::vector<G4Track*> trackVector;

        for (int i = 0; i < trajectoryContainer->size(); i++) {
            G4Trajectory* trajectory = (*trajectoryContainer)[i];
            if (trajectory->GetParentID() == 0) {
                trackVector.push_back(trajectory->GetTrack());
            }
        }

        // Generate scintillation photons for each track
        for (int i = 0; i < trackVector.size(); i++) {
            G4Track* track = trackVector[i];
            if (track->GetDefinition() == G4Electron::ElectronDefinition() && track->GetKineticEnergy() > 1.0 * keV) {
                GenerateScintillationPhotons(track);
            }
        }
    }

private:
    void GenerateScintillationPhotons(G4Track* aTrack)
    {
        // Generate the scintillation photons
        G4StepPoint* preStepPoint = aTrack->GetStep()->GetPreStepPoint();
        G4ThreeVector pos = preStepPoint->GetPosition();
        G4ThreeVector dir = preStepPoint->GetMomentumDirection();

        G4int nScintillationPhotons = fScintProcess->GetNumPhotons();
        G4double scintillationTime = fScintProcess->GetMeanLifeTime();
        G4double scintillationRiseTime = fScintProcess->GetScintillationRiseTime();
        G4double scintillationDecayTime = fScintProcess->GetScintillationDecayTime();

        for (int i = 0; i < nScintillationPhotons; i++) {
            // Generate a random direction for the photon
            G4double cosTheta = 2.0 * G4UniformRand() - 1.0;
            G4double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
            G4double phi = 2.0 * M_PI * G4UniformRand();
            G4ThreeVector photonDir(sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta);

            // Generate a random time for the photon
            G4double photonTime = preStepPoint->GetGlobalTime() + G4RandExponential::shoot(scintillationTime);

            // Generate a random energy for the photon
            G4double photonEnergy = G4RandGauss::shoot(1.0 * eV, 0.1 * eV);

            // Create the photon and set its properties
            G4DynamicParticle* photon = new G4DynamicParticle(fOpticalPhoton, photonDir, photonEnergy);
            photon->SetParentID(aTrack->GetTrackID());

            // Add the photon to the event
            G4PrimaryVertex* vertex = new G4PrimaryVertex(pos, photonTime);
            vertex->SetPrimary(new G4PrimaryParticle(photon));
            fEvent->AddPrimaryVertex(vertex);

            // Add the photon to the waveform
            fScintSD->AddWaveformValue(photonTime, photonEnergy, scintillationRiseTime, scintillationDecayTime);
        }
    }

    G4Event* fEvent;
    G4ParticleDefinition* fOpticalPhoton;
    G4Scintillation* fScintProcess;
    G4OpBoundaryProcess* fOpBoundaryProcess;
    ScintillatorSD* fScintSD;
};

其中，`ScintillatorSD` 类是一个用户定义的类，用于在Geant4模拟中保存闪烁体上的光信号波形。下面是这个类的代码：

class ScintillatorSD : public G4VSensitiveDetector
{
public:
    ScintillatorSD(const G4String& name) : G4VSensitiveDetector(name), fWaveformSize(0)
    {
        for (int i = 0; i < WAVEFORM_SIZE; i++) {
            fWaveform[i] = 0.0;
        }
    }

    virtual ~ScintillatorSD() {}

    virtual void Initialize(G4HCofThisEvent*)
    {
        fWaveformSize = 0;
        for (int i = 0; i < WAVEFORM_SIZE; i++) {
            fWaveform[i] = 0.0;
        }
    }

    virtual G4bool ProcessHits(G4Step* aStep, G4TouchableHistory*)
    {
        return true;
    }

    void AddWaveformValue(G4double time, G4double amplitude, G4double riseTime, G4double decayTime)
    {
        // Calculate the index of the first waveform bin that overlaps with this photon
        G4int startIndex = std::max(0.0, std::floor((time - 5.0 * decayTime) / WAVEFORM_TIME_STEP));

        // Calculate the index of the last waveform bin that overlaps with this photon
        G4int endIndex = std::min(WAVEFORM_SIZE - 1.0, std::ceil((time + 5.0 * riseTime) / WAVEFORM_TIME_STEP));

        // Add the waveform values for this photon to the appropriate bins
        for (int i = startIndex; i <= endIndex; i++) {
            G4double t = i * WAVEFORM_TIME_STEP;
            fWaveform[i] += amplitude * std::exp(-(t - time) / decayTime) * (1.0 - std::exp(-(t - time) / riseTime));
        }

        // Update the waveform size
        fWaveformSize = std::max(fWaveformSize, endIndex + 1);
    }

    void InitializeWaveform()
    {
        fWaveformSize = 0;
        for (int i = 0; i < WAVEFORM_SIZE; i++) {
            fWaveform[i] = 0.0;
        }
    }

    G4int GetWaveformSize() const
    {
        return fWaveformSize;
    }

    G4double GetWaveformValue(G4int index) const
    {
        return fWaveform[index];
    }

private:
    static const G4int WAVEFORM_SIZE = 10000;
    static const G4double WAVEFORM_TIME_STEP = 0.01 * ns;

    G4double fWaveform[WAVEFORM_SIZE];
    G4int fWaveformSize;
};

在这个示例代码中，我们使用了一个简单的敏感探测器 `ScintillatorSD`，它保存了一个大小为10000的数组，用于存储光信号波形。每个光子在产生时都会调用 `AddWaveformValue` 函数，该函数将光子的波形值添加到数组中。在模拟结束时，我们将数组中的所有值写入一个名为 "waveform.txt" 的文本文件中，以便后续分析。