1

我正在尝试使用 Pol01 示例来检查 G4_Si 散射的初始偏振光子康普顿的角分布。
我必须在 SteppingAction.cc 中包含几行代码，这些代码将打印出诸如康普顿散射角和光子的最终动能之类的信息。

模拟似乎适用于水平和垂直偏振光子、线性偏振 +- 45 度和右圆偏振 (RCP) 光子。但是，当我尝试对左圆偏振 (LCP) 光子运行模拟时，我在将信息打印到终端窗口的过程中遇到了分段错误。

在宏代码中，这是我发现的

### RCP (works)
# /gun/polarization 0. 0. 1.

### LCP ( Segmentation fault )
/gun/polarization 0. 0. -1.

### vertical ( works )
# /gun/polarization -1. 0. 0.

### horizontal ( works )
# /gun/polarization 1. 0. 0.

### LP ( 45 degrees ) (works)
# /gun/polarization 0. 1. 0.

### LP ( -45 degrees ) ( works )
# /gun/polarization 0. -1. 0.

/gun/particle gamma
#
/gun/energy 511 keV

但是，在此代码上方有以下说明

/polarization/volume/set theBox 0. 0. 0.08

如果我设置0.08为0.，那么我不会得到分段错误。我想帮助弄清楚为什么只有 LCP 光子会导致分段错误？

仅供参考，我包含在 Pol01 中的唯一额外代码（在 SteppingAction.cc 中）如下

  // Inside this method you should access the track from the step object:
  G4Track *aTrack = aStep->GetTrack();

  // Then you should get the  KineticEnergy of the track if the process is 
  // "Compton scattering" and if the volume remains the same in the next step. 
  // See the code lines below:
  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();

  const G4VProcess* CurrentProcess = preStepPoint->GetProcessDefinedStep();
  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 == "pol-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 = 13 (angular spread (in radians) as a function of scattering angle)
        // analysisManager->FillH1(13, thetaRad);

        // 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 z = pos.z();   

        // Check if volume is sensitive detector
        G4String fromVolume = aTrack->GetNextVolume()->GetName();
        // see https://geant4-forum.web.cern.ch/t/get-particle-info-in-steppingaction-if-particle-pass-through-specific-volume/2863/2
        // for possible solution for finding which volume the gamma is in
        // G4String SDname = aStep->GetPostStepPoint()->GetSensitiveDetector()->GetName();
        // FOR CHECKING
        std::cout << "fromVolume: "                 << fromVolume
                  << "    Particle: "               << ParticleName
                  << "    Process: "                << StepProcessName
                  << "    Energy: "                 << EkinUnits 
                  << "    Angle: "                  << thetaRad/degree << " deg." << std::endl;
      }
    }
  }

感谢您的时间。彼得

4

0 回答 0