Scintillation light

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

in a post before (Scintillation light - #14 by faca87) you mentioned the issue that it did not print anything, as if n_proc was zero. that is weird now that it is always 10…
do you strictly copy paste the code from your program to your posts?

Hello @weller, it didn’t print anything because for error, I wrote the printing line in other “if”…

I fixed that issue and this is my stepping action but I still don’t get scintillation photons
B1SteppingAction.cc (48.0 KB)

and with

do you mean that counting does not work, or that they are not tracked (not present in the simulation at all)?

what were the results of testing this code snippet:

Hello, @weller

I wrote the code (I used the filter “if” to be sure that I’m in the scintillator)

if(NextVol && ThisVol->GetName()=="PlasticScintillator" && NextVol->GetName()=="PlasticScintillator") { 
    const G4ParticleDefinition* particleDef = step->GetTrack()->GetDynamicParticle()->GetParticleDefinition();
    	if(particleDef == opticalphoton) G4cerr << "I am a photon :)" << G4endl;

and here the log
log.txt (2.8 KB)

as you can see, the sentence I am a photon :slight_smile: has been written only 3 times.

Moreover, the code (this code is out of the filter)

s

tatic 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;
      		}

doesn’t find scintillation photons

Here the stepping action

B1SteppingAction.cc (48.2 KB)

then to me it seems the best way for you is to directly track the photons, and don’t use the second way.

the number of photons can depend on so many aspects, you will have to investigate :slight_smile:

Hello, then, what should I do?

maybe simplify your code to the bare minimum, make it work there, slowly add more „complex“ things and check if things still work as they should.
you could also use the linked example as a start, from the filenames it looks like you took B1 instead of OpNovice2!?

Hello @weller

Yes, I started from B1 example! This is because I’m implementing the scintillation to my simulation…I’m building the simulation step -by-step and my supervisor said me to start from B1

Searching on the forum, I found this @dsawkey 's message Track of optical photon - #5 by dsawkey he said that in the macro of the OpNovice2 there are some command lines needed to use the scintillation I guess the lines are:

/process/optical/processActivation Cerenkov      false
/process/optical/processActivation OpAbsorption  true
/process/optical/processActivation OpBoundary    false
/process/optical/processActivation Scintillation true
/process/optical/processActivation OpRayleigh    false
/process/optical/processActivation OpMieHG       false
/process/optical/processActivation OpWLS         false

/process/optical/scintillation/setEnhancedTimeConstants true
/process/optical/scintillation/setYieldFactor 10
#/process/optical/scintillation/setExcitationRatio .5
/process/optical/scintillation/setByParticleType true
/process/optical/scintillation/setTrackInfo      false
/process/optical/scintillation/setFiniteRiseTime false
/process/optical/scintillation/setStackPhotons   true

/process/optical/scintillation/setTrackSecondariesFirst true

but if I add these lines I also don’t get the sentence
if(particleDef == opticalphoton) G4cerr << "I am a photon :)" << G4endl;

the creation of photons is a statistical process, depending on type and energy of particles inside the scintillator volume, and also (very important) the geometry of the scintillator and the scintillating properties that you defined.
how many particles reach the scintillator when you don’t get any photons?

this is why I suggested to make the simulation as simple as possible, to be able to really understand what is happening.
to make things work in this case you really only need 1 volume with your scintillating material, and nothing else.
make the primary particles start right in front of the scintillator volume, point source, pencil beam, so all primaries make it into the scintillator.

for my personal taste, the stepping action should be empty except for the things you want to test…

Hello @weller I simplified the simulation with only the scintillator directly hitted by the beam,

This is the reduced simulation WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free

but I still don’t get scintillation photons

played around with the code a little bit, especially with the stuff from OpNovice2 example, and must admit, I could not yet find out why this code in the stepping action does not produce any photon counts, even though photons are created and being stacked and traced:

    // this demonstrates use of GetNumPhotons()
    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();
      }
    }

probably only a single line of configuration missing / different than in the example!? Did you build the example and does it work in your installed geant4 version?

I found some presentation slides where a similar code snippet is slightly different, you could give it a try. Don’t know if that is supposed to work fine with the current version of geant4?

G4int photons = 0;
G4SteppingManager* fpSteppingManager = 
  G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager();
G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
if (stepStatus != fAtRestDoItProc) {
  G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector(); 
  size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops(); 
  for (size_t i3 = 0; i3 < MAXofPostStepLoops; i3++) {
    if ((*procPost)[i3]->GetProcessName() == "Scintillation") {
      G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3];
      photons += proc1->GetNumPhotons();
    }
  }
}

edit: not helpful, unfortunately :frowning:

What version of geant4 are you using? Just switched the test simulation to version 11, and voila:

G4WT0 > In this step, 1 scintillation photons were produced.
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)
G4WT0 > In this step, 1 scintillation photons were produced.
G4WT0 > I am a photon :)
G4WT0 > In this step, 1 scintillation photons were produced.
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)
G4WT0 > In this step, 1 scintillation photons were produced.
G4WT0 > I am a photon :)
G4WT0 > I am a photon :)

this is with the code block from my last post:

    G4int photons = 0;
    G4SteppingManager* fpSteppingManager =
      G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager();
    G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
    if (stepStatus != fAtRestDoItProc) {
      G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector();
      size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops();
      for (size_t i3 = 0; i3 < MAXofPostStepLoops; i3++) {
        if ((*procPost)[i3]->GetProcessName() == "Scintillation") {
          G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3];
          photons += proc1->GetNumPhotons();
        }
      }
    }
    if(photons > 0)
    {
        G4cerr << "In this step, " << photons
               << " scintillation photons were produced." << G4endl;
    }

remember now to look up the correct scintillation properties for your material :wink:

Hello @weller, thank you. I’m having problems with Geant4 11 …I will open a topic for it and then I will try your code!

EDIT : @weller did to you use my full code …or did you use your simulation? because I’m having problems to use my simulation by Geant4 11 because it can’t find g4analysis.hh as I wrote here Geant4 11 on cvmf

Hello @weller I finally can run on G4 11
but I can’t get scintillation photons using this code in the stepping action


	G4int photons = 0;
				G4SteppingManager* fpSteppingManager = 
  				G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager();
				G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
				if (stepStatus != fAtRestDoItProc) {
  					G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector(); 
  					size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops(); 
  					for (size_t i3 = 0; i3 < MAXofPostStepLoops; i3++) {
    					if ((*procPost)[i3]->GetProcessName() == "Scintillation") {
      						G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3];
      						photons += proc1->GetNumPhotons();
    					}
  					}	
				}	
      			if(photons > 0){ 
			  		G4cout << "In this step, "<< photons << " scintillation photons were produced." << G4endl;
			  			analysisManager->FillNtupleDColumn(8,0, photons);
					       		analysisManager->AddNtupleRow(8);  
      			}

B1SteppingAction.cc (47.9 KB)

this in the detector construction

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

B1DetectorConstruction.cc (30.0 KB)

and this in the running macro


/process/optical/processActivation Cerenkov      false
/process/optical/processActivation OpAbsorption  true
/process/optical/processActivation OpBoundary    false
/process/optical/processActivation Scintillation true
/process/optical/processActivation OpRayleigh    false
/process/optical/processActivation OpMieHG       false
/process/optical/processActivation OpWLS         false

run1.txt (7.3 KB)

This is the log file

log.txt (108.7 KB)

It was a very minimal simulation, that was in principle just the mandatory classes plus the stepping action. and there is only 1 G4Box in the whole world :wink:

I feel this file is extremely hard to debug, especially from outside… what are all these if(someFlag) doing there? how can you be sure that there is no logical mistake? how can the forum community verify that this is the case? → make it simple! copy the whole application to a new folder, and strip out everything that is not needed to debug the issue.

my bet now would be on: fPlastScintFlag == false, hence ScintMat = world_mat and there all these G4MaterialProperties are not defined.

on top: how many particles are passing through the scintillator? what is the thickness the particles effectively see? if it is only few particles and a thin thickness → maybe you need many more primaries to be lucky and see a photon…