0

我有一个简单的几何结构，由一块薄硅板（我的敏感探测器）组成，我在其中发射一束 511 keV 光子的准直光束。我只想选择那些用敏感探测器散射过一次康普顿的光子。

我试图在我的 SensitiveDetector.cc（称为 SensitiveDetectorCPET03.cc）中放置一个条件来执行此操作，但没有成功。我希望使用 G4Step 类来告诉我主光子是否是康普顿散射的。

我的代码部分（在 SensitiveDetectorCPET03.cc 中）尝试应用此条件如下

（HitCPET04.hh 是 myHit 脚本）

G4bool SensitiveDetectorCPET04::ProcessHits(G4Step* aStep, G4TouchableHistory*){

HitCPET04* newHit = new HitCPET04();

我想访问 GetProcessDefinedStep() 因为我的印象是这会告诉我主光子和敏感探测器（如果有的话）之间发生了什么物理过程。所以我首先介绍这行代码。

G4StepPoint* postStepPoint = aStep->GetPostStepPoint();

现在我尝试创建一个名为“process”的常量指针，我希望它指向
存储在 GetProcessDefinedStep() 中的信息，希望它可以告诉我是什么过程导致我的光子散射。所以我写了这行代码

```const G4VProcess *process = postStepPoint->GetProcessDefinedStep();```

如果我可以写一个条件，比如..

if( *process == "compt" ){

        fHitsCollection->insert( newHit );
        // get analysis manager```
        auto analysisManager = G4AnalysisManager::Instance();

        G4double En = KineticEnergy/keV;
        G4double theta = acos( 2- 511.0/En );        
        // fill ntuple
        analysisManager->FillNtupleDColumn(0, KineticEnergy/keV); 
        analysisManager->FillNtupleDColumn(1, theta/deg); 

        analysisManager->AddNtupleRow(0); 
        std::cout << postStepPoint << std::endl;
    }       
return true;

}

那么我想我会很好，但这不起作用。我的编译器对我大喊大叫

没有运算符 "==" 匹配这些操作数 -- 操作数类型是： const G4VProcess == const char [6]

我不知道如何使用此错误消息来修复我的代码。

我是 Geant4 的新手，所以我提前道歉，我的知识非常有限。我想如何设置条件以仅记录那些用敏感探测器散射一次康普顿的光子。

感谢您花时间阅读我的请求。

此致

彼得

4

1 回答 1

0

在这个出色的解释的帮助下，我已经解决了这个问题。 https://www.researchgate.net/post/How_to_stop_particles_tracking_in_GEANT4

下面的直方图显示了当我将角度分布与 Klein-Nishina 理论进行比较时得到的一致性。（标题应为 0.5 毫米立方硅。）

[克莱因-仁科理论与模拟][1]

（以下是对我的 ProccessHits 代码的调整+一些解释。我希望这对其他人有帮助。）

G4bool SensitiveDetectorCPET04::ProcessHits(G4Step* aStep, G4TouchableHistory*) { // 实例化 G4Analysismanager

auto analysisManager = G4AnalysisManager::Instance();

// cf https://www.researchgate.net/post/How_to_stop_particles_tracking_in_GEANT4 // 在此方法中，您应该从 step 对象访问轨道：

G4Track *aTrack = aStep->GetTrack();

// 如果过程是 // “康普顿散射”并且音量在下一步中保持不变，那么您应该获得轨道的 KineticEnergy。// 请参见下面的代码行：

G4StepPoint* preStepPoint = aStep->GetPreStepPoint();


if (CurrentProcess != 0)
 {

// To implement the next line, one needs to include " #include "G4VProcess.hh" "
// (See SteppingAction.cc in TestEm6)
    
const G4String & StepProcessName = CurrentProcess->GetProcessName();
// Get track ID for current processes

    G4double TrackID = aStep->GetTrack()->GetTrackID();
       if( (StepProcessName == "compt")&&(TrackID==1) ) 
       {

    // processing hit when entering the volume
    G4double kineticEnergy = aStep->GetTrack()->GetKineticEnergy();

    auto EkinUnits = G4BestUnit(kineticEnergy, "Energy");

    // Get energy in units of keV              
    G4double En = kineticEnergy/keV;

    // Determine angle in radians
   G4double thetaRad = std::acos( 2- 511/En );

    // Place condition to reject Compton scattering at 0 rad.
    // (This is to avoid nan results when dividing through by 2pisin(theta))

       if ( (thetaRad > 0) && (thetaRad < pi))
      {
      // fill histogram id = 0 (angular spread (in radians) as a function of scattering angle)        
      analysisManager->FillH1(0, thetaRad); 
// (This produces the data in the above graph) ---> The theoretical curve comes from my root script.

      // fill of histogram id = 1 (angular spread (in radians) as a function of solid angle)
      // c.f. TestEm5 (TrackingAction.cc) and also (corrections at)
      // https://geant4-forum.web.cern.ch/t/physics-validation/2893

      G4double dthetaRad1  = ( analysisManager->GetH1Width(1) );

      G4double unit1   = analysisManager->GetH1Unit(1);

      G4double weight1 = (unit1)/( twopi*std::sin(thetaRad)*dthetaRad1 );

      analysisManager->FillH1(1, thetaRad, weight1); **GRAPH NOT SHOWM**
   
      // get particle name that's undergone the Compton scattering
    
      G4String ParticleName = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
      
// Find position (c.f. /examples/extended/medical/GammaTherapy) in CheckVolumeSD.cc

      G4ThreeVector pos = aStep->GetPreStepPoint()->GetPosition();

      // G4double x = pos.x();
      // G4double y = pos.y();
      G4double z = pos.z();   

      // Fill NTuples 
      analysisManager->FillNtupleDColumn(0, En/keV); 
      analysisManager->FillNtupleDColumn(1, thetaRad/deg); 
      // analysisManager->FillNtupleDColumn(2, x/um);
      // analysisManager->FillNtupleDColumn(3, y/um);
      analysisManager->FillNtupleDColumn(2, z/um);
      analysisManager->AddNtupleRow(0); 
      // FOR CHECKING
      std::cout << "fromVolume: "                 << fromVolume
                // << "    unit: "                   << unit
                << "    Particle: "               << ParticleName
                << "    Process: "                << StepProcessName
                << "    Energy: "                 << EkinUnits 
                << "    Angle: "                  << thetaRad/degree << " deg." << std::endl;
      //          << "    Scattering site: "        << pos  << std::endl; 
      }
    }
  }
return true;
}

Best
Peter


[1]: https://i.stack.imgur.com/WsKu1.jpg
于 2020-09-29T12:59:54.083 回答