Failed To Make Something That would Work Like An SIPM in GEANT4

Geometry, Fields and Transportation
1

Hello dear GEANT users,

I have a scintillator created in geant4 I also want to make an SIPM for that scintillator or at least something that would give the same results as it does, but I was not succesfull at make an such SIPM like system and I need advice. I share the parts of my .cc files relevant to SIPM creation. But I am told that they are hard to read so let me explain them. First of all I wrapped the scintillator with aluminum except for 6mm x 6mm part of it(SIPM will be put to that part). Then I made a 6mm x 6mmx0.41mm silicon cube to resemble an Sensl 30065. Then I created a logical border surface with silicon and scintillator like this:

//------------------------------Optical Surface of sipm-scint------------------ 
    
    G4OpticalSurface* OpSurface_for_sipm = new G4OpticalSurface("scintillator_sipm_Surface");
    
    //G4LogicalBorderSurface* Surface = new G4LogicalBorderSurface("Scint1Surface", phys_scint1, phys_surface_scint1, OpSurface);
    
    OpSurface_for_sipm->SetType(dielectric_dielectric);
    //OpSurface->SetFinish(polishedfrontpainted);
    OpSurface_for_sipm->SetFinish(polished);
    OpSurface_for_sipm->SetModel(unified);
    OpSurface_for_sipm->SetSigmaAlpha(0.1*deg); // no known
    
    std::vector<G4double> pp_of_Si           = {  2.0 * eV, 3.47 * eV };//taken from: /examples/extended/optical/wls/src/WLSMaterials.cc
    std::vector<G4double> rindex_of_Si       = {1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46, 1.46 }; 
    //std::vector<G4double> reflectivity = { 1.0, 1.0 };
    //std::vector<G4double> efficiency   = { 1.0, 1.0 };

    G4MaterialPropertiesTable* opSurfaceProperty_for_sipm = new G4MaterialPropertiesTable();

    opSurfaceProperty->AddProperty("REFLECTIVITY", photonEnergy, efficiency);
    opSurfaceProperty->AddProperty("EFFICIENCY", photonEnergy, reflectivity);
    opSurfaceProperty_for_sipm->AddProperty("RINDEX", photonEnergy, rindex_of_Si);
    //opSurfaceProperty_for_sipm->AddProperty("ABSLENGTH", pp_of_Si, absClad);
    OpSurface_for_sipm->SetMaterialPropertiesTable(opSurfaceProperty_for_sipm);
    
    G4LogicalBorderSurface* Surface_sipm = new G4LogicalBorderSurface("Scint1Surface", phys_scint1, phys_sipm1, OpSurface_for_sipm);

//------------------------------End of Optical Surface of Simp-scint-------------------------

then I kind of tried to replicate what an SIPM does and few extra stuff with SIPMSD.cc and SIPMHit.cc . At the end of SIPMSD.cc I tried to print out all of scintillators results but nothing gets printed even the numbers I put for control. This is how I print stuff out(I tried root files before this that did not work either):

G4bool SipmSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
{
  // energy deposit
  
    G4double h = 4.13566733e-15*eV*s;
    G4double c = 299792458.*m/s;
    
    G4cout << "1" << G4endl;
  
  if(aStep->GetTrack()->GetDefinition()->GetParticleName() != (G4String)"opticalphoton") return false; //idea stolen from: LXePMTSD.cc
  
  G4cout << "2" << G4endl;
  
  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
  
  G4cout << "3" << G4endl;
  
  G4TouchableHistory* touchable = (G4TouchableHistory*)(preStepPoint->GetTouchable());
  
  G4cout << "4" << G4endl;
  
  G4VPhysicalVolume* pv = aStep->GetPostStepPoint()->GetPhysicalVolume();
  
  if(pv->GetName() != "PVsipm1")return false;
  
  G4cout << "5" << G4endl;
  
  G4double edep = aStep->GetTotalEnergyDeposit();
  
  G4double totalEn = aStep->GetTrack()->GetTotalEnergy();
  
  G4double wavelength = (h*c)/totalEn;

  G4double preTime = aStep->GetPreStepPoint()->GetGlobalTime();

  G4double postTime = aStep->GetPostStepPoint()->GetGlobalTime();
  
  G4cout << "6" << G4endl;
  
  G4cout << "wavelength of photon: " << wavelength << "\n Post time of photon: " << postTime << "\n Pre time of photon: " << preTime << G4endl;  
  
  G4cout << "7" << G4endl;
  
  G4cout.flush();
  //G4cout << 1 << G4endl;

  //if (edep==0.) return false;

  auto newHit = new SipmHit();

  G4cout << "8" << G4endl;

  newHit->SetTrackID  (aStep->GetTrack()->GetTrackID());
  
  G4cout << "9" << G4endl;
  
  G4int number=-1;
  G4int copyno=-1;
  if(collectionName[0] == "DetCollection1") number=0;
  //if(collectionName[0] == "DetCollection2") number=1;
  //if(collectionName[0] == "DetCollection3") number=2;
  
  //copyno = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber();
  
  G4cout << "10" << G4endl;
  
  newHit->SetChamberNb(number);
  newHit->SetEdep(edep);
  //newHit->SetCopyNumber(copyno);
  newHit->SetPos (aStep->GetPostStepPoint()->GetPosition());
  newHit->SetTime(postTime);
  newHit->SetTime(preTime);
  newHit->SetWavelength(wavelength);
  newHit->SetTotalEnergy(totalEn);

  G4cout << "11" << G4endl;

  fHitsCollection->insert( newHit );

  G4cout << "12" << G4endl;
  
  aStep->GetTrack()->SetTrackStatus(fStopAndKill);

  G4cout << "13" << G4endl;

  return true;
}

Well at this point I am out of ideas and require assistance.

DetectorConstruction.cc (15.2 KB)
SipmSD.cc (6.1 KB)

_Geant4 Version:_11.1
_Operating System:_Ubuntu 20.04
_Compiler/Version:_9.4.0
_CMake Version:_3.25.1

2

It appears that I can not share all my code at once since I am a new user so this is the rest of the code:
SipmHit.cc (3.0 KB)
PhysicsList.cc (8.1 KB)

3

Are you calling ProcessHists during your stepping action?

If you are never getting here:

G4cout << “1” << G4endl;

Then your ProcessHits function is never called. I would assume you want to call this during your stepping action, either during every step (I wouldn’t) or after some checking, such as if your PV is your SiPM

4

Hello lodyms, first of all thanks for the answer,

No I do not call ProcessHits during stepping action and I dont know how to do that. my stepping action is like this:

/// \file B1/src/SteppingAction.cc
/// \brief Implementation of the B1::SteppingAction class

#include “SteppingAction.hh”
#include “G4SteppingManager.hh”
#include “EventAction.hh”
#include “DetectorConstruction.hh”
#include “G4AnalysisManager.hh”
#include “RunAction.hh”
#include “G4Step.hh”
#include “G4Event.hh”
#include “G4RunManager.hh”
#include “G4LogicalVolume.hh”
#include “G4Electron.hh”
#include “G4Positron.hh”
#include <G4AntiNeutrinoMu.hh>

#include “G4Types.hh”
#include “G4UserStackingAction.hh”

namespace B1
{

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

SteppingAction::SteppingAction(EventAction* eventAction)
: fEventAction(eventAction)
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void SteppingAction::UserSteppingAction(const G4Step* step)
{
  if (!fScoringVolume) {
    const auto detConstruction = static_cast<const DetectorConstruction*>(
      G4RunManager::GetRunManager()->GetUserDetectorConstruction());
    fScoringVolume = detConstruction->GetScoringVolume();
  }

  G4AnalysisManager* man = G4AnalysisManager::Instance();
  // get volume of the current step
  G4String volume
    = step->GetPreStepPoint()->GetTouchableHandle()
      ->GetVolume()->GetLogicalVolume()->GetName();

  G4double toten = step->GetPreStepPoint()->GetTotalEnergy();
    
  G4int trackid = step->GetTrack()->GetTrackID();
  
  G4int parentid = step->GetTrack()->GetParentID();
  
  G4double time = step->GetTrack()->GetGlobalTime()/CLHEP::ns;
  
  //G4String name = step ->GetTrack()->GetParticleDefinition()->GetParticleName();
  
  //G4cout << name << G4endl;
  
/*
  if(trackid == 1)
  {
    man->SetNtupleMerging( 1 );
    man->FillNtupleDColumn(0,0,toten);
    man->AddNtupleRow(); 
  }*/
  
    // Retrieve the secondary particles
    
    // for more info /examples/advanced/radioprotection
    
    G4SteppingManager*  steppingManager = fpSteppingManager;
    G4TrackVector * fSecondary = steppingManager -> GetfSecondary();

  /*	
   for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
     { 
       // Retrieve the info about the generation of secondary particles 
       G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); // volume where secondary was generated 
       
       G4String secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();  // name of the secondary
       
       G4double secondaryParticleKineticEnergy =  (*fSecondary)[lp1] -> GetKineticEnergy(); // kinetic energy
       
       G4String process = (*fSecondary)[lp1]-> GetCreatorProcess()-> GetProcessName(); //what kind of physics happened, stolen from: ./advanced/brachytherapy/src/BrachySteppingAction.cc

       G4cout << "\n secondary's name : " << secondaryParticleName << "\n secondary id: " << trackid << " \n secondaries time: "<< time << " \n secondary's parent id: " << parentid << "\n Secondary's Logical Volume: " << volume << "\n Process Name: " << process <<G4endl; 
    /*if (step -> GetTrack() -> GetTrackID() == 1 && step->GetTrack()->GetDefinition() == G4Electron::Electron() ) {
  
  G4cout << "e-" <<G4endl;
  
  }*/

    
  // check if we are in scoring volume
  //if (volume != fScoringVolume) return;

  // collect energy deposited in this step
  G4double edepStep = step->GetTotalEnergyDeposit();
  fEventAction->AddEdep(edepStep);


}
  

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

What should I write/include in SteppingAction to get ProcessHits’ printings working. Also yes in ProcessHits it is state to check if PV is SIPM.

5

Hello again,

After your reply I messed around the examples to find something that does what I want to do, I came across this:

     for ( G4int icount=0; icount<nofHits; icount++ ){ (*fHitsCollection)[icount]->Print();
     G4cout << (*fHitsCollection)[icount] << G4endl;}

I also created an EndofEvent function to use what is stated above. Did not change anything unfortunately.

6

So you have all of this code in your function ProcessHits in your namespace SipmSD. But you need to call it during a step.

So you need to call SipmSD::ProcessHits(step,G4TouchableHistory*) during your stepping action, or the code that is in SiPMSD will never be executed.

Also, you could just have all of that processing in your Stepping action.

In the end, what you want to do is for each step, check if it is an optical photon hitting your SiPM. If it is, then you need to record that in some histogram.

7

Hello again and thanks for the answer,

I will do what you just said in stepping action right now. Also I do not know how to call SipmSD::ProcessHits on SteppingAction how can it be done surely including headher files and writing SipmSD::ProcessHits(); is not enough.

8

You attach a sensitive detector (SD) instance to a volume. The Geant4 tracking code will automatically call yourSD::ProcessHits() at each step in that volume, provided you attached it correctly and wrote the class correctly.

You don’t generally want to call an SD from your own SteppingAction, since you’re duplicating effort (and double counting if you did register the SD correctly).

1 Like
9

Hello and thanks for the reply,

I must say I am confused on what I should do to see my results, returning to the main point of the post, I would like to hear recommendations from you as well.

10

Mkelsey is correct, you do not want to call SiPMSD::ProcessHits() directly (although… you could write a custom tracking/detector system this way).

I have not worked with sensitive detectors so I did not recognize that is what you are trying to do. You need to link your SiPM SD class to your SiPM volume so that when a particle steps through it, the processhits function is called (mkelsey can correct me if that is wrong)

11

Hello and thanks for the answer,

I am doing linking stuff you mentioned with this, in the detectorconstruction.cc code:
(at least that is what I think I am doing.)

void DetectorConstruction::ConstructSDandField()
{

  // Sensitive detectors

  G4String SipmChamberSD1name = "/SipmSD1";
  auto aSipmSD1 = new SipmSD(SipmChamberSD1name, "DetCollection1");
  G4SDManager::GetSDMpointer()->AddNewDetector(aSipmSD1);
  SetSensitiveDetector("logic_sipm1", aSipmSD1, true);
  
}
12

Can you paste in your SetSensitiveDetector() function? It looks like you have that function set up to take the name of the LV, rather than an explicit pointer to a (saved) LV object.

13

I don’t get what you mean by paste SetSensitiveDetector() function. It is not a function I write by myself or copied. It is something I find in example B2. In example B2 it used like this and there were no extra stuff neither in .cc or .hh of detector construction(as far as I know);

void DetectorConstruction::ConstructSDandField()
{
  // Sensitive detectors

  G4String trackerChamberSDname = "B2/TrackerChamberSD";
  auto aTrackerSD = new TrackerSD(trackerChamberSDname, "TrackerHitsCollection");
  G4SDManager::GetSDMpointer()->AddNewDetector(aTrackerSD);
  SetSensitiveDetector( fLogicChamber,  aTrackerSD );
14

I see! I’ve learned something new today. That SetSensitiveDetector() function is provided by the base G4VUserDetectorConstruction class for convenience (which means I didn’t need to create buffers to store the LVs to be made sensitive).

So then the question comes back to whether your SD is written to record the things you want it to record. If you want to save the information to histograms or N-tuples, you can do that with the G4Analysis tools. You could instatiate ROOT class objects to store information, and do the file handling yourself in your code. You could use simple cout statements and save the information to a dumb text file.

1 Like
15

It looks like he implements his output, information storing, etc in the function

SipmSD::ProcessHits()

I assume the process is to set up your volume as an SD so that whenever a track/step passes through this volume, ProcessHits() is called. It looks like this is done in this line:

auto aSipmSD1 = new SipmSD(SipmChamberSD1name, “DetCollection1”);
G4SDManager::GetSDMpointer()->AddNewDetector(aSipmSD1);
SetSensitiveDetector(“logic_sipm1”, aSipmSD1, true);

1 Like
16

I indeed intended to do that, yet I have no idea how to solve my problems.

17

You have these print statements in your ProcessHits function. Are you saying these do not run? Are you sure that you have particles that are taking steps in your volume?

18

Hello once again,

I say I have no output also I see particles and their interaction with scintillator and scintillation light from gui.

19

Do you have your entire code in a git-hub repository?

20

I have it in a drive:https://drive.google.com/drive/folders/1RdKZXaNexfmtCPAS5vuwHUpJYU8wmpca?usp=drive_link

If a drive does not suit you of course I can share it on github.

Thanks for reply.