How select hits composed of Compton scattering events?

Particles, Track, Event, Run and Biasing
1

Hi.
I am really new at Geant4, so I apologize in advance that my knowledge is very limited, but I am trying very hard to understand how to use Geant4.

I have a simple geometry consisting of a thin silicon plate (my sensitive detector) in which I fire a collimated beam of 511 keV photons at it. I would like to select just the Compton scattered photons that have scattered once inside the sensitive detector.

I attempted to place a condition in my SensitiveDetector.cc (called SensitiveDetectorCPET03.cc) which would do this, but without success. I was hoping to use the aStep class to tell me if the primary photon was Compton scattered.

Here is the section of my code which I am attempting to apply this condition

G4bool SensitiveDetectorCPET03::ProcessHits(G4Step* aStep, G4TouchableHistory*)
{
G4double TrackID = aStep->GetTrack()->GetTrackID();
G4String ParticleName = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
G4double KineticEnergy = aStep->GetPostStepPoint()->GetKineticEnergy();
G4ThreeVector pos = aStep->GetPostStepPoint()->GetPosition();

if (KineticEnergy == 0.511) return false;
    
HitCPET03* newHit = new HitCPET03();
   
if( KineticEnergy < 0.511 && TrackID == 2)
// if( KineticEnergy < 0.511)
    {
        fHitsCollection->insert( newHit );
        // get analysis manager
        auto analysisManager = G4AnalysisManager::Instance();

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

        analysisManager->AddNtupleRow(0); 
    } 
          
return true;

}

The graph I get when I plot counts verses scatering angle does not match the Klein-Nishina formula.

I thought I found a away to include a condition to select just compton scattering events in a set of slides I found here --> https://conferences.fnal.gov/g4tutorial/g4cd/Slides/IEEE/readout.pdf

In slide 12 it shows a line of code which reads
G4StepPoint* preStepPoint = aStep->GetPreStepPoint();

I thought that perhaps I could use this G4SetPoint to determine if the photon that collided with the detector was due to Compton scattering by writing something like this

G4StepPoint* process = aStep->GetPreStepPoint()->GetProcessDefinedStep();
then write something like
if (process == “compt”)

However, if I include
“G4StepPoint* process = aStep->GetPreStepPoint()->GetProcessDefinedStep();”
into my SensitiveDetectorCPET03.cc script my Visual Studio (running in Ubuntu virtual box) doesn’t like it and screams at me with the message
“a value of type “const G4VProcess *” cannot be used to initialize an entity of type “G4StepPoint *””

I do not know what this means.

I don’t know whether this may be of interest but at the top of my SensitiveDetectorCPET03.cc file I have the following header files included.

#include “AnalysisCPET03.hh”
#include “SensitiveDetectorCPET03.hh”
#include “Randomize.hh”
#include “G4SDManager.hh”
#include “G4RandomDirection.hh”
#include “G4SystemOfUnits.hh”

I would feel very grateful if someone out there could point me in the right direction in terms of finding a way to record just Compton scattering hits.

Thank you for your consideration.
Best regards
Peter

1 Like
2

This is exactly what you want, but you need to fix your code. The compiler message is telling you exactly what to fix. You’ve tried to assign the output of the GetProcessDefinedStep() function to a variable of type G4StepPoint*. If you look at G4StepPoint.hh, you’ll discover (as the compiler error tells you) that that function returns a pointer of type const G4VProcess*. Your variable needs to be of that type.

3

Dear mkelsey,
Thank you for taking time out to assist with my coding issue, and confirming that I am on the right track by using the code “G4StepPoint* process = aStep->GetPreStepPoint()->GetProcessDefinedStep();”
which simplifies things a bit for me. I do also realize the compiler is telling me exactly what to fix, but I do not understand what it means.

I have looked in the G4StepPoint.hh header file, as you suggested, and on line 161 I found this line of code
const G4VProcess* GetProcessDefinedStep() const;

Now, keep in mind I am only half way through my C++ training (learning about functions and function overloading) and I am about 2 chapters away before I am ready to learn about Objects and Classes in greater detail.
My best guess is that the this part; aStep->GetPostStepPoint()->GetProcessDefinedStep() is a constant string member. So I tried making the pointer postStepPoint to this member a constant pointer such that

const G4StepPoint *postStepPoint = aStep->GetPostStepPoint()->GetProcessDefinedStep();

But this is wrong because the compiler is still yelling at me with the message
a value of type “const G4VProcess *” cannot be used to initialize an entity of type “const G4StepPoint *”

I do not have enough c++ experience to understand how to use this message to fix my coding issue.

So I am still stuck. But thanks for letting me know that I am at least on the right track.
Best regards
Peter

4

C++ is a “strongly typed” language (unlike Python, for example). If you are assigning the output of a function to a variable, that variable must be declared with the same type as what the function returns. In C++ code, the type declaration appears before the name of the function.

Thus,

double someFunction() { return 7.6; }

. . .

std::string myValue = someFunction();

will not compile, because someFunction() is returning a double but you’re trying to put that return value into a std::string. You are in the same situation here with GetProcessDefinedStep().

5

I have solved the problem, with the help of this excellent explanation.
https://www.researchgate.net/post/How_to_stop_particles_tracking_in_GEANT4

The histogram below shows the agreement I get when I compare the angular distribution compared to Klein-Nishina theory. The title should read 0.5 mm cube of Si.
(The adjustment to my ProccessHits code + some explanation follows. I hope this can help others.)

KNvsSIM_0pt5mmx5mmSi_SingleScatter
KNvsSIM_0pt5mmx5mmSi_SingleScatter998Ă—912 29.8 KB

G4bool SensitiveDetectorCPET04::ProcessHits(G4Step* aStep, G4TouchableHistory*)
{
  // Instantiate G4Analysismanager
  auto analysisManager = G4AnalysisManager::Instance();

  // c.f. https://www.researchgate.net/post/How_to_stop_particles_tracking_in_GEANT4
  // 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 == "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) ---> I implement the theoretical curve in my root script.

        // fill of histogram id = 1 (angular spread (in radians) as a function of soild 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