Inconsistent Particle Hits on Sensitive Detector

[kmbartel]

_Geant4 Version:_v11.3.0
_Operating System:_macOS Sequoia 15.2
_Compiler/Version:_v16.0.0
_CMake Version:_v3.31.2

---

Hello! I’ve been using Geant4 to model particle production from a 120 GeV proton beam incident on an iron beam dump, and have been attempting to record particle name and 4-momentum for any particles passing through 2 detectors at different places in the simulation geometry. I’ve recently noticed something a little odd that I was wondering if anyone had any clarifications on: the number of particles that appear to pass through each detector don’t seem to line up with the output on the visualizer, especially for the first detector. Here’s one example run below, with the visualizer and the corresponding detector hits, where there are significantly fewer particles being recorded than there should be according to the visualizer:

Below is the code that I use for both detectors. My aim was to have each detector record the particle name and 4-momentum of the particle when it first hit the detector and both print it out and save it to a .csv file, though it doesn’t seem as though that is the case based on the above image. What could be causing this discrepancy and how can I fix it? Any help would be greatly appreciated!

#include "detector.hh"
#include "G4ParticleDefinition.hh"
#include "G4AnalysisManager.hh"


MySensitiveDetector::MySensitiveDetector(G4String name) : G4VSensitiveDetector(name)
{}

MySensitiveDetector::~MySensitiveDetector()
{}

G4bool MySensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *ROHist)
{
    //new KE, mom. stuff from G4 forums
    //check to make sure track isn't leaving the worl vol. first (returns physical vol. as null ptr)
    G4bool worldboundary = aStep->GetPostStepPoint()->GetStepStatus()==fWorldBoundary;
    if (!worldboundary){

        G4LogicalVolume* post_volume = aStep->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
        G4String vol_name = post_volume->GetName();
        G4bool boundary = aStep->GetPreStepPoint()->GetStepStatus()==fGeomBoundary;

        if ((vol_name == "Detector") && (boundary)){
            //get kinetic energy, momenta, particle name
            G4double KE = aStep->GetPostStepPoint()->GetKineticEnergy();
            G4double p_x = aStep->GetPostStepPoint()->GetMomentum().x(); 
            G4double p_y = aStep->GetPostStepPoint()->GetMomentum().y();
            G4double p_z = aStep->GetPostStepPoint()->GetMomentum().z();
        
            G4String particleName = aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetParticleName();

            //print out the 4-mom.
            G4cout << "Particle Name, 4-mom.: " << vol_name << ", " << particleName << ", " << "(" << KE << ", " << p_x << ", " << p_y << ", " << p_z << ")" << G4endl;

            //record the 4-mom.
            G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
            analysisManager->FillNtupleSColumn(0, 0, vol_name);
            analysisManager->FillNtupleSColumn(0, 1, particleName);
            analysisManager->FillNtupleDColumn(0, 2, KE);
            analysisManager->FillNtupleDColumn(0, 3, p_x);
            analysisManager->FillNtupleDColumn(0, 4, p_y);
            analysisManager->FillNtupleDColumn(0, 5, p_z);
            analysisManager->AddNtupleRow(0);
        }
        
    }

     return 1;

}

#include "seconddetector.hh"
#include "G4ParticleDefinition.hh"
#include "G4AnalysisManager.hh"


MySecondSensitiveDetector::MySecondSensitiveDetector(G4String name) : G4VSensitiveDetector(name)
{}

MySecondSensitiveDetector::~MySecondSensitiveDetector()
{}

G4bool MySecondSensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *ROHist)
{
    //new KE, mom. stuff from G4 forums
    //check to make sure track isn't leaving the worl vol. first (returns physical vol. as null ptr)
    G4bool worldboundary = aStep->GetPostStepPoint()->GetStepStatus()==fWorldBoundary;
    if (!worldboundary){

        G4LogicalVolume* post_volume = aStep->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
        G4String vol_name = post_volume->GetName();
        G4bool boundary = aStep->GetPreStepPoint()->GetStepStatus()==fGeomBoundary;
        
        if ((vol_name == "DetectorTwo") && (boundary)){
            //get kinetic energy, momenta, particle name
            G4double KE = aStep->GetPostStepPoint()->GetKineticEnergy();
            G4double p_x = aStep->GetPostStepPoint()->GetMomentum().x(); 
            G4double p_y = aStep->GetPostStepPoint()->GetMomentum().y();
            G4double p_z = aStep->GetPostStepPoint()->GetMomentum().z();
        
            G4String particleName = aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetParticleName();

            //print out the 4-mom.
            G4cout << "Particle Name, 4-mom.: " << vol_name << ", " << particleName << ", " << "(" << KE << ", " << p_x << ", " << p_y << ", " << p_z << ")" << G4endl;

            //record the 4-mom.
            G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
            analysisManager->FillNtupleSColumn(0, 0, vol_name);
            analysisManager->FillNtupleSColumn(0, 1, particleName);
            analysisManager->FillNtupleDColumn(0, 2, KE);
            analysisManager->FillNtupleDColumn(0, 3, p_x);
            analysisManager->FillNtupleDColumn(0, 4, p_y);
            analysisManager->FillNtupleDColumn(0, 5, p_z);
            analysisManager->AddNtupleRow(0);
        }
        
    }

     return 1;

}

[pico]

I would to take one step back and try to see the bigger picture here:

It looks to me like you are are assigning the *SensitiveDetector* to the world volume? And you are actually doing it twice(?). One where you check for "Detector" and another for
"DetectorTwo". Is that correct?

I would probably assign the same MySensitiveDetector to the two logical volumes for the detectors and join the logic into one single class.

I’m also a simple person… I would shoot a geantino and see that it CONSISTENTLY crosses my two detectors, otherwise there might be a flaw in your logic

/Pico

[kmbartel]

Hi–thanks so much for helping out! To answer your first question, I define two different sensitive detectors (detector.cc and second detector.cc) and incorporate them into the overall program by defining two logical volumes ("Detector" and "DetectorTwo") in my detector construction. Here’s the code for the detector:

#include "construction.hh"

MyDetectorConstruction::MyDetectorConstruction()
{}

MyDetectorConstruction::~MyDetectorConstruction()
{}

G4VPhysicalVolume *MyDetectorConstruction::Construct()
{
    G4NistManager *nist = G4NistManager::Instance();
    G4Material *worldMat = nist->FindOrBuildMaterial("G4_AIR");
    
    //volumes are always half-volumes (i.e. divide the actual volume by two)
    //old world volume:
    G4Box *solidWorld = new G4Box("solidWorld", 2.5*m, 3.5*m, 13.*m);

    //G4Box *solidWorld = new G4Box("solidWorld", 2.5*m, 3.5*m, 6.*m);
    
    G4LogicalVolume *logicWorld = new G4LogicalVolume(solidWorld, worldMat, "logicWorld");
    
    G4VPhysicalVolume *physWorld = new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), logicWorld, "physWorld", 0, false, 0, true);
    
    //we need 7 shapes for the most basic detector sketch:
    
    //SHAPE 1: FMag and Beam Dump (iron)
    G4Material *beamdumpMat = nist->FindOrBuildMaterial("G4_Fe");
    //old position:
    G4ThreeVector pos1 = G4ThreeVector(0, 0, -9.5*m);

    //G4ThreeVector pos1 = G4ThreeVector(0, 0, -2.5*m);
    
    //new dimensions (1 meter depth):
    G4Box *solidShape1 = new G4Box("BeamDump", 2.3*m, 2.5*m, 0.5*m);

    //original dimensions:
    //G4Box *solidShape1 = new G4Box("BeamDump", 2.3*m, 2.5*m, 2.5*m);
    
    //G4LogicalVolume *logicShape1 = new G4LogicalVolume(solidShape1, beamdumpMat, "BeamDump");
    logicShape1 = new G4LogicalVolume(solidShape1, beamdumpMat, "BeamDump");
    
    new G4PVPlacement(nullptr,  // no rotation
      pos1,                     // at position
      logicShape1,              // its logical volume
      "BeamDump",               // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //set the color and style
    G4VisAttributes *beamdumpVisAtt = new G4VisAttributes(G4Colour(1., 1., 0.));
    beamdumpVisAtt->SetForceSolid(true);
    logicShape1->SetVisAttributes(beamdumpVisAtt);
    
    //DETECTOR
    G4Material *detectorMat = nist->FindOrBuildMaterial("G4_AIR");
    //old detector placement
    G4ThreeVector posDetect = G4ThreeVector(0, 0, -6.625*m);

    //G4ThreeVector posDetect = G4ThreeVector(0, 0, -0.625*m);

    //G4ThreeVector posDetect = G4ThreeVector(0, 0, -9.625*m);
    
    G4Box *solidDetector = new G4Box("Detector", 2.5*m, 3.5*m, 0.125*m);
    
    //if you don't want a segfault, make sure to not use the * operator here
    //G4LogicalVolume *logicDetector = new G4LogicalVolume(solidDetector, detectorMat, "Detector");
    logicDetector = new G4LogicalVolume(solidDetector, detectorMat, "Detector");
    
    new G4PVPlacement(nullptr,  // no rotation
      posDetect,                // at position
      logicDetector,            // its logical volume
      "Detector",               // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //SHAPE 2: Station 1 (silicon)
    G4Material *stat1Mat = nist->FindOrBuildMaterial("G4_Si");
    G4ThreeVector pos2 = G4ThreeVector(0, 0, -6.*m);
    
    G4Box *solidShape2 = new G4Box("Station1", 1.*m, 1.*m, 0.5*m);
    
    G4LogicalVolume *logicShape2 = new G4LogicalVolume(solidShape2, stat1Mat, "Station1");
    
    new G4PVPlacement(nullptr,  // no rotation
      pos2,                     // at position
      logicShape2,              // its logical volume
      "Station1",               // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //set the color and style
    G4VisAttributes *stat1VisAtt = new G4VisAttributes(G4Colour(1., 0., 1.));
    stat1VisAtt->SetForceSolid(true);
    logicShape2->SetVisAttributes(stat1VisAtt);
    
    //SHAPE 3: KMag (iron)
    G4Material *kmagMat = nist->FindOrBuildMaterial("G4_Fe");
    G4ThreeVector pos3 = G4ThreeVector(0, 0, -2.5*m);
    
    G4Box *solidShape3 = new G4Box("initKMag", 2.*m, 2.*m, 1.5*m);
    
    G4Box *holeBox = new G4Box("Hole", 1.*m, 1.*m, 1.5*m);
    
    G4SubtractionSolid *finalsolidShape3 = new G4SubtractionSolid("solidKMag", solidShape3, holeBox);
    
    G4LogicalVolume *logicShape3 = new G4LogicalVolume(finalsolidShape3, kmagMat, "KMag");
    
    new G4PVPlacement(nullptr,  // no rotation
      pos3,                     // at position
      logicShape3,              // its logical volume
      "KMag",                   // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //set the color and style
    G4VisAttributes *kmagVisAtt = new G4VisAttributes(G4Colour(0., 1., 1.));
    kmagVisAtt->SetForceSolid(true);
    logicShape3->SetVisAttributes(kmagVisAtt);
    
    //SHAPE 4: Station 2 (silicon)
    G4Material *stat2Mat = nist->FindOrBuildMaterial("G4_Si");
    G4ThreeVector pos4 = G4ThreeVector(0, 0, 0.5*m);
    
    G4Box *solidShape4 = new G4Box("solidStation2", 1.5*m, 1.5*m, 0.5*m);
    
    G4LogicalVolume *logicShape4 = new G4LogicalVolume(solidShape4, stat2Mat, "Station2");
    
    new G4PVPlacement(nullptr,  // no rotation
      pos4,                     // at position
      logicShape4,              // its logical volume
      "Station2",               // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //set the color and style
    G4VisAttributes *stat2VisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5));
    stat2VisAtt->SetForceSolid(true);
    logicShape4->SetVisAttributes(stat2VisAtt);
    
    //SHAPE 5: Station 3 (silicon)
    G4Material *stat3Mat = nist->FindOrBuildMaterial("G4_Si");
    G4ThreeVector pos5 = G4ThreeVector(0, 0, 4.5*m);
    
    G4Box *solidShape5 = new G4Box("solidStation3", 1.5*m, 1.5*m, 0.5*m);
    
    G4LogicalVolume *logicShape5 = new G4LogicalVolume(solidShape5, stat3Mat, "Station3");
    
    new G4PVPlacement(nullptr,  // no rotation
      pos5,                     // at position
      logicShape5,              // its logical volume
      "Station3",               // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //set the color and style
    G4VisAttributes *stat3VisAtt = new G4VisAttributes(G4Colour(1., 0., 0.));
    stat3VisAtt->SetForceSolid(true);
    logicShape5->SetVisAttributes(stat3VisAtt);

    //DETECTOR 2
    G4Material *detectorTwoMat = nist->FindOrBuildMaterial("G4_AIR");
    
    G4ThreeVector posDetectTwo = G4ThreeVector(0, 0, 5.625*m);
    
    G4Box *solidDetectorTwo = new G4Box("DetectorTwo", 2.5*m, 3.5*m, 0.125*m);
    
    //if you don't want a segfault, make sure to not use the * operator here
    //G4LogicalVolume *logicDetector = new G4LogicalVolume(solidDetector, detectorMat, "Detector");
    logicDetectorTwo = new G4LogicalVolume(solidDetectorTwo, detectorTwoMat, "DetectorTwo");
    
    new G4PVPlacement(nullptr,  // no rotation
      posDetectTwo,             // at position
      logicDetectorTwo,         // its logical volume
      "DetectorTwo",            // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //SHAPE 6: Absorber Wall (lead)
    G4Material *absorberMat = nist->FindOrBuildMaterial("G4_Pb");
    G4ThreeVector pos6 = G4ThreeVector(0, 0, 7.25*m);
    
    G4Box *solidShape6 = new G4Box("solidAbsorber", 2.*m, 2.*m, 0.5*m);
    
    G4LogicalVolume *logicShape6 = new G4LogicalVolume(solidShape6, absorberMat, "Absorber");
    
    new G4PVPlacement(nullptr,  // no rotation
      pos6,                     // at position
      logicShape6,              // its logical volume
      "Absorber",               // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //set the color and style
    G4VisAttributes *absorbVisAtt = new G4VisAttributes(G4Colour(0., 0., 1.));
    absorbVisAtt->SetForceSolid(true);
    logicShape6->SetVisAttributes(absorbVisAtt);
    
    //SHAPE 7: Station 4 (silicon)
    G4Material *stat4Mat = nist->FindOrBuildMaterial("G4_Si");
    G4ThreeVector pos7 = G4ThreeVector(0, 0, 9.5*m);
    
    G4Box *solidShape7 = new G4Box("solidStation4", 2.*m, 2.*m, 1.5*m);
    
    G4LogicalVolume *logicShape7 = new G4LogicalVolume(solidShape7, stat4Mat, "Station4");
    
    new G4PVPlacement(nullptr,  // no rotation
      pos7,                     // at position
      logicShape7,              // its logical volume
      "Station4",               // its name
      logicWorld,               // its mother  volume
      false,                    // no boolean operation
      0,                        // copy number
      true);                    // overlaps checking
    
    //set the color and style
    G4VisAttributes *stat4VisAtt = new G4VisAttributes(G4Colour(0., 1., 0.));
    stat4VisAtt->SetForceSolid(true);
    logicShape7->SetVisAttributes(stat4VisAtt);
    
    return physWorld;
}

void MyDetectorConstruction::ConstructSDandField()
{
    MySensitiveDetector *sensDet = new MySensitiveDetector("SensitiveDetector");
    
    logicDetector->SetSensitiveDetector(sensDet);

    MySecondSensitiveDetector *sensDetTwo = new MySecondSensitiveDetector("SecondSensitiveDetector");
    
    logicDetectorTwo->SetSensitiveDetector(sensDetTwo);
    
    G4MagneticField *magField = new G4UniformMagField(G4ThreeVector(1.8*tesla, 0., 0.));
    G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
    fieldMgr->SetDetectorField(magField);
    fieldMgr->CreateChordFinder(magField);

    G4bool allLocal = true;
    logicShape1->SetFieldManager(fieldMgr, allLocal); //here insert the name of the logical vol. with the B-field
}

Do you think that my assigning two different sensitive detectors this way could be causing the issues that I’m seeing?

I also tried your suggestion of using a geantino, and the visualizer showed that the geantino track passed through both detectors for the several events I ran, so it doesn’t seem as though there are any issues with particle tracks.

Thanks so much for your help!

[jrellin]

This looks like you are trying to just count flux through the two surfaces. Sensitive Detector might not be what you want and may cause unexpected behavior since it will only “trigger” if there is an interaction and/or energy deposited in an air volume in the first step in that volume. I am not surprised that you aren’t seeing many counts.

I would just pass or grab the logical volumes and check that is where you are. Then the particle just needs a step, any step, into the volume for your checks to work.

[kmbartel]

Thank you so much for the advice! I’m a little confused on how to go about recording particle data without a Sensitive Detector. Could you give me a bit more guidance on how to go about that, or point me to an example that does something similar?

[jrellin]

You basically have it. Move the code from your Sensitive Detector to your stepping action instead. You no longer need Sensitive Detector files for these two air volumes. All that will change is:

# Now inside Stepping Action
if ((vol_name == "Detector") && (boundary))

to

# Now inside Stepping Action
if (((vol_name == "Detector") || (vol_name == "Detector2")) && (boundary))

Then in your NTuple I would just add a field that is the Detector “ID” that can be 1 or 2 depending on which volume you are in. Then everything is cleaner and you should see a lot more particles.

[kmbartel]

This was exactly correct–I’m now seeing the appropriate number of particles. Thanks so much!

[jrellin]

No problem. One last thing you might want to add for your check is to check the sign of the dot product of the particles momentum with the normal vector of the air “volumes”. From your image this would be the vector pointing to the “right” in the first image you shared. Some particles, especially neutrons, can backscatter which would cause you to double count them since they’d be coming first from the right and then from the left of an air volume after colliding with some other detector solid. Checking the dot product is not a complete fix for this but will work most of the time.