Scintillation light

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

I added it in my main file

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

exampleB1.cc (20.5 KB)

it should be my

G4Material* scintillator = nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
 G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
   MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
   scintillator->SetMaterialPropertiesTable(MPT);

in the detector construction
B1DetectorConstruction.cc (29.7 KB)

I added those lines in the stepping action
B1SteppingAction.cc (47.1 KB)

When I run, I get many Process: Scintillation printed by the line
G4cerr << "Process: " << aStep->GetPostStepPoint()->GetProcessDefinedStep()-

but I don’t get any scintillation photons (that should be no phyisically possible…if there is scintillation light…there are photons)…

if it can help, this is the full simulation code WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free

ok, think I got it:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html?#example-physproc-5

the material presumably needs some more properties being defined.

    auto* silicium = nistManager->FindOrBuildMaterial("G4_Si");

    std::vector<G4double> energy     = {2.034*eV, 3.*eV, 4.136*eV};
    std::vector<G4double> rindex     = {1.3435, 1.351, 1.3608};
    std::vector<G4double> absorption = {344.8*cm, 850.*cm, 1450.0*cm};

    G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();

    // property independent of energy
    MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
    // properties that depend on energy
    MPT->AddProperty("RINDEX", energy, rindex);
    MPT->AddProperty("ABSLENGTH", energy, absorption);

    silicium->SetMaterialPropertiesTable(MPT);

in the stepping action:

static G4ParticleDefinition* opticalphoton =
            G4OpticalPhoton::OpticalPhotonDefinition();

    const G4ParticleDefinition* particleDef =
            aStep->GetTrack()->GetDynamicParticle()->GetParticleDefinition();

    if(particleDef == opticalphoton)
        G4cerr << "I am a photon :)" << G4endl;

...
G4WT0 > Process: Scintillation
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)
...

edit: in a way this makes sense… without a defined refractive index, what should the wavelength be? defaulting to zero would yield division by zero.

Hello @weller, thank you

I tried to add the code

  1. DetectorConstruction
 //Scintillation process
   std::vector<G4double> energy     = {2.034*eV, 3.*eV, 4.136*eV};
    std::vector<G4double> rindex     = {1.3435, 1.351, 1.3608};
    std::vector<G4double> absorption = {344.8*cm, 850.*cm, 1450.0*cm};

   G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
   MPT->AddProperty("RINDEX", energy, rindex);
    MPT->AddProperty("ABSLENGTH", energy, absorption);
    
   MPT->AddConstProperty("SCINTILLATIONYIELD", 11800./MeV);
   scintillator->SetMaterialPropertiesTable(MPT);

and steppingaction

> static G4ParticleDefinition* opticalphoton = G4OpticalPhoton::OpticalPhotonDefinition();
> 				const G4ParticleDefinition* particleDef = track->GetDynamicParticle()->GetParticleDefinition();
>  			auto proc_man = track->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();
>         			analysisManager->FillNtupleDColumn(8,0, n_scint);
>         			analysisManager->AddNtupleRow(8);  
> 					}
>     			}
>       		if( n_scint > 0){ 
> 			  G4cout << "In this step, "<< n_scint << " scintillation photons were produced." << G4endl;
>       			}

but I still don’t get scintillation photons

B1SteppingAction.cc (47.2 KB)
B1DetectorConstruction.cc (30.0 KB)

what if you switch back to the variant without iterating over the proc_vec?
and what is values is n_proc when there is scintillation?

in which way?

and what is values is n_proc when there is scintillation?
[/quote]
I tried to print
G4cout << "n_proc==" <<n_proc<< G4endl;
both with and without the filter

if(step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "Scintillation"){

and it always writes