How select hits composed entirely of Compton scattering events?

Question

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 want to select just those photons that have Compton scattered once with 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 G4Step class to tell me if the primary photon was Compton scattered.

The section of my code (in SensitiveDetectorCPET03.cc) which attempts to apply this condition is as follows

(HitCPET04.hh is myHit script)

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

HitCPET04* newHit = new HitCPET04();

I want to access GetProcessDefinedStep() because I am under the impression that this will tell me what physical process happened between my primary photon and sensitive detector (if any). So I start by introducing this line of code.

G4StepPoint* postStepPoint = aStep->GetPostStepPoint();

Now I attempt to create a const pointer called "process" which I hope points to the
information stored in GetProcessDefinedStep() hoping that it can tell me what process caused my photon to scatter. So I write this line of code

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

If I could then write a condition like..

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;

}

then I think I would be good to go but this does not work. My compiler screams at me with the message

no operator "==" matches these operands -- operand types are: const G4VProcess == const char [6]

I do not know how to use this error message to fix my code.

I am really new at Geant4, so I apologize in advance that my knowledge is very limited. I would like to how to set the condition to record only those photons that have one time Compton scattered with the sensitive detector.

Thank you for taking the time to read my request.

Best regards

Peter

Answer 1

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.)

[Klein-Nishina theory vs Simulation][1]

(The adjustment to my ProccessHits code + some explanation follows. I hope this can help others.)

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();

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