Scintillation light

Hello, I simulated particles entering a scintillator and I get their features (number, energy etc) and the total deposited energy in the scintillator.

Is it possibile also simulating how many of them can produce scintillation light in the scintillator?

There’s nothing specific to scintillation to do this. I suggest writing a UserSteppingAction and checking if a particle creates scintillation photons. Examples are in the extended/optical examples, e.g. here:
https://geant4.kek.jp/lxr/source/examples/extended/optical/OpNovice2/src/SteppingAction.cc#L324

Hello @dsawkey thank you…I tried to follow the example…but surely I did a mistake.
Firstly I want to explain you my goal.

In november 2021 it has been measured the radiation hitting a scintillator connected to an ArduSiPM. The beam (1GeV electron beam) entered and exit a steel vacuum chamber by two mylar chamber; In front of a third vacuum chamber forming an angle of 22.5 degree with the beam axis, there was the scintillator to measure the radiation outgoing the vacuum chamber. The counting rate measured by was 4.410^5 count/s.
I simulated the setup, and I estimated the counting rate of particles ENTERING the scintillator…it is 1.5
10^7 count/s…
This is the geant4 pictorial view (green=vacuum chamber, blue=mylar windows, grey=scintillator)

Obiviously in the 1.510^7 particles there are many of them that ArduSiPM can’t measure. For example, to be detected by the SiPM, photons entering in the scintillator must produce elctron-positron pairs and the latter must produce scintillation light. If I try to apply the E>1.2MeV to photons entering the scintillator (i.e. the threshold to produce electron-positron pair), the number of photons is reduced from 1.210^7 to 7.1*10^5)…
therefore, now I want to look by the simulation how many of the particles entering the scintillator produce scintillation light in it (i.e. can in principle be detected by the SiPM)

I wrote in the stepping action

if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="PlasticScintillator") { 
						int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
			 			analysisManager->FillNtupleDColumn(6,0, pdg);
			 			double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
			 			analysisManager->FillNtupleDColumn(6,1, kinEnergy);
			 			double MomDirx = step->GetTrack()->GetMomentumDirection().x();
			 			analysisManager->FillNtupleDColumn(6,2, MomDirx);
			 			double MomDiry = step->GetTrack()->GetMomentumDirection().y();
			 			analysisManager->FillNtupleDColumn(6,3, MomDiry);
			 			double MomDirz = step->GetTrack()->GetMomentumDirection().z();
			 			analysisManager->FillNtupleDColumn(6,4, MomDirz);
			 			double Vertx = step->GetTrack()->GetVertexPosition().x();
			 			analysisManager->FillNtupleDColumn(6,5, Vertx);
						double Verty = step->GetTrack()->GetVertexPosition().y();
						analysisManager->FillNtupleDColumn(6,6, Verty);
						double Vertz = step->GetTrack()->GetVertexPosition().z();
						analysisManager->FillNtupleDColumn(6,7, Vertz);
						double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
						analysisManager->FillNtupleDColumn(6,8, Mom);
						double Dirxsp = step->GetTrack()->GetPosition().x();
			 			analysisManager->FillNtupleDColumn(6,9, Dirxsp); 
			 			double Dirysp = step->GetTrack()->GetPosition().y();
			 			analysisManager->FillNtupleDColumn(6,10, Dirysp); 
			 			double Dirzsp = step->GetTrack()->GetPosition().z();
			 			analysisManager->FillNtupleDColumn(6,11, Dirzsp); 
			 			double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
						analysisManager->FillNtupleDColumn(6,12, Momx);
						double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
						analysisManager->FillNtupleDColumn(6,13, Momy);
						double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
						analysisManager->FillNtupleDColumn(6,14, Momz);
						double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(6,15, xprime);
						double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
						analysisManager->FillNtupleDColumn(6,16, yprime);
						auto proc_man = step->GetTrack()->GetParticleDefinition()->GetProcessManager();
						G4String proc_name = step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
    					G4int n_scint = 0;
						 if(proc_name == "Scintillation"){
						 	G4Scintillation* scint = new G4Scintillation("Scintillation");
						 	n_scint    = scint->GetNumPhotons();
					       	analysisManager->FillNtupleDColumn(6,17, n_scint);
				            }
			   			analysisManager->AddNtupleRow(6);  
					}

but you see…I get 0 scintillation photons (it was just to try, so I simulated only 10^6 primary electrons)

B1SteppingAction.cc (46.0 KB)

Probably you’re not generating any opticalphotons. Are you sure you’ve got all the material properties set correctly? Run with /tracking/verbose 2 and grep for opticalphoton.

Hello @dsawkey

Now simulated 30*10^6 primary electrons, but I still don’t get scintillation photons…then I think the problem is the code …

In my DetecorConstruction

G4Material* scintillator = nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
	G4Box* PlastScint;
G4Material* ScintMat = scintillator;
PlastScint = new G4Box("PlastScint", 0.5*fScint_sizeX, 0.5*fScint_sizeY, 0.5*fScint_sizeZ);
     logicScint = new G4LogicalVolume(PlastScint, ScintMat, "PlasticScintillator");

Then you see…I used the pre-set G4_PLASTIC_SC_VINYLTOLUENE

B1DetectorConstruction.cc (29.4 KB)

I tried to enable
/tracking/verbose 2
but in few seconds the file on which Geant4 writes the tracking was of approximately 1GB…so I don’t think that it’s possible to enable it…or is there a way to write the tracking only for opticalphoton???

EDIT

Hello @dsawkey
I tried to print the process in this way

if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="PlasticScintillator") { 
	G4String proc_name = step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
						 G4cout << "Process " << proc_name <<  G4endl;
    					G4int n_scint = 0;
						 if(proc_name == "Scintillation"){
						 	G4Scintillation* scint = new G4Scintillation("Scintillation");
						 	n_scint    = scint->GetNumPhotons();
					       	analysisManager->FillNtupleDColumn(6,17, n_scint);
				            }

I simulated 300*10^6 primary electrons, but you see…I don’t get scintillation process
log.txt (360.3 KB)

maybe is it because I’m looking the step from the world to the scintillator?

But I also tried to look for steps into the scintillator i.e.

if(NextVol && ThisVol->GetName()=="PlasticScintillator" && NextVol->GetName()=="PlasticScintillator") {

but in this way it looks like that after few time, the simulation doesn’t continue the run (i.e. there aren’t error message, but it doesn’t process new events)…sorry if I’m in a hurry, but I need to make this work before than january 29

just to be sure: this is with multiple runs, right? per run you can only simulate about 2e9 events.

Sorry @weller! it’s a typo! I simulated 300*10^6 events

this is all in the bracket of the if-case

if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="PlasticScintillator") {

when the track crosses a geometric boundary, I don’t think the GetProcessDefinedStep will ever be scintillation, it will always be the one for crossing boundaries.

Thank you @weller, then what should I do to compare results of the simulation (i.e. 1.210^7 particles entering the scintillator) with the measurement by the ArduSiPM (i.e. 4.410^5 particles)?
What should I do to get scintillation?

As I wrote, I tried to look only in the scintillator i.e.

if(NextVol && ThisVol->GetName()=="PlastiScintillator" && NextVol->GetName()=="PlasticScintillator") {
	 G4cout << "Process " << proc_name <<  G4endl;
    					G4int n_scint = 0;
						 if(proc_name == "Scintillation"){
						 	G4Scintillation* scint = new G4Scintillation("Scintillation");
						 	n_scint    = scint->GetNumPhotons();
					       	analysisManager->FillNtupleDColumn(6,17, n_scint);
				            }
}

but after some runs it doesn’t make anything…i.e. there aren’t errors messages, but it looks like that the simulation doens’t simulate all the primary particles

Here the example…it simulates 1.2*10^6 events, then stops without errors
log.txt (3.0 KB)

is this a typo? missing a c there…

the log looks ok to me, note that the different workers are processing the events independently. so after the 1.2e6 event start the remaining work load is simply finished.

fyi: Just checked because I was unsure: the process for crossing geometric boundaries is indeed CoupledTransportation.

Yes…I checked the stepping action, there is the “c”

but for other run, I read alll the current runs and the log ends with the root file closing!

Yest it is the crossing geometry CoupledTransportation

Now I tried to write

	auto proc_man = step->GetTrack()->GetParticleDefinition()->GetProcessManager();
						G4String proc_name = step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();
					//	 G4cout << "Process " << proc_name <<  G4endl;
    					G4int n_scint = 0;
						 if(proc_name == "Scintillation"){
						 	G4Scintillation* scint = new G4Scintillation("Scintillation");
						 	n_scint    = scint->GetNumPhotons();
					       	analysisManager->FillNtupleDColumn(8,0, n_scint);
					       		analysisManager->AddNtupleRow(8);  	
				            }

out of the “if”… and to save the scintillation photons in a new root ttree…I simulated 3*10^6 electrons, but I don’t get scintillation photons…
I really don’t understand…

have you compared your code to the OpNovice2 example that was linked here before? mabe you could try to adapt your code to how it was implemented there?

the main difference I see at first glance is that they iterate over a vector

G4ProcessVector* proc_vec = proc_man->GetPostStepProcessVector(typeDoIt);

instead of using

step->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName();

don’t know if that could solve your issue?

yes @weller! I just did it! I wrote in the stepping action

auto proc_man = step->GetTrack()->GetDynamicParticle()->GetParticleDefinition()->GetProcessManager();
    		G4ProcessVector* proc_vec = proc_man->GetPostStepProcessVector(typeDoIt);
    		G4int n_proc = proc_vec->entries();
			G4int n_scint = 0;
    		for(G4int i = 0; i < n_proc; ++i) {
      			G4String proc_name = (*proc_vec)[i]->GetProcessName();
			    if(proc_name.compare("Scintillation") == 0) {
        			auto scint = (G4Scintillation*) (*proc_vec)[i];
			        n_scint    = scint->GetNumPhotons();
      				}
      				G4cout << "In this step, " << n_scint << " scintillation photons were produced." << G4endl;
    			}

Then it would write something about n_scint…i.e. I didn’t write if n_scint>0,…so even if n_scint==0 it should write
In this step, 0 scintillation photons were produced.

but look the log…it didn’t write anything

log.txt (5.0 KB)

B1SteppingAction.cc (46.8 KB)

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html#pre-defined-properties
“Starting in version 11.0, some optical material properties are defined in Geant4.”

→ scintillation seems to not be defined in most / all the predefined materials. I think you have to add this yoursel :slight_smile:

example: Physics Processes — Book For Application Developers 11.0 documentation

Hello @wellr,

  1. firstly I fixed an issue…it didn’t write the sentence, because I wrote the G4couet in other if …sorry! now I fixed it

B1SteppingAction.cc (46.9 KB)

  1. I added the scintillation property in the detector construction
  //Scintillation process
   G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
   MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
   scintillator->SetMaterialPropertiesTable(MPT);

B1DetectorConstruction.cc (29.7 KB)

I simulated 3*10^6 primary electrons, but it didn’t get scintillator photons

even if 39 particles entered in the scintillator

PS
This line
G4int n_proc = proc_vec->entries();
produces this warinng

C:\B1\src\B1SteppingAction.cc(784,41): warning C4267: 'inizializzazione': conversione da 'size_t' a 'G4int'. Possibile
perdita di dati. [C:\B1\B1-build\exampleB1.vcxproj]

PPS. to have more particles, I set the scintillator material to the target…then 3*10^6 primary electrons (800MeV) cross a scintillator material…but I still don’t get scintillator photons…

what physics list do you use in your code?
does it somehow resemble this from the examples:

 G4VModularPhysicsList* physicsList = new FTFP_BERT;
 physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
 G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
  
 physicsList->RegisterPhysics(opticalPhysics);
 runManager->SetUserInitialization(physicsList);
 

the warning about size_t type not being G4int can be ignored.

This is the physicslist:
PhysicsList.cc (8.4 KB)

Can you upload the file to see where adding the code? in the example
examples\extended\optical\OpNovice2\src

I don’t see the physiclist file…

it is the few lines that I quoted for you. You find it here:
https://geant4.kek.jp/lxr/source/examples/extended/optical/OpNovice2/OpNovice2.cc#L59

no trace of G4OpticalPhysics in this file. so that could explain why there is no optical photons :slight_smile:

Thank you…then I modified in the main file
from

 // Physics list
  G4VModularPhysicsList* physicsList = new QBBC;
  physicsList->SetVerboseLevel(1);
  runManager->SetUserInitialization(physicsList);

to

 G4VModularPhysicsList* physicsList = new FTFP_BERT;   
  physicsList->ReplacePhysics(new G4EmStandardPhysics_option4());
  G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
  	physicsList->RegisterPhysics(opticalPhysics);
  physicsList->SetVerboseLevel(1);
  runManager->SetUserInitialization(physicsList);

exampleB1.cc (20.4 KB)

but I still don’t get scintillation photon, even trying to set the target with scintillator material…

that is indeed weird. I just implemented the steps as described, and it seems to work.

just for completeness, the steps:

  1. Register G4OpticalPhysics to the physics list by adding this line:
physicsList->RegisterPhysics(new G4OpticalPhysics());
  1. Add Scintillation properties to the material in question:
    auto* silicium = nistManager->FindOrBuildMaterial("G4_Si");
    G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
    MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
    silicium->SetMaterialPropertiesTable(MPT);
  1. Print/Filter out Scintillation events in the stepping action:
    if(aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "Scintillation")
        G4cerr << "Process: " << aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() << G4endl;
  1. Rebuild code :wink:
  2. Check output:
...
G4WT1 > Process: Scintillation
...